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1 Introduction 

1.1 Summary 


Gravity gradiometers in satellites have received a lot of attention because of their ability to 
measure a combination of the local gravity gradient, plus instrument rotation effects. After 
a series of measures to isolate the gradient, a global mesh of measurements can be combined 
to determine the planetary external gravity potential, a result of great importance in 
geophysics and geodynamics. 

In 1993, Mr. Seymor Kant, the sponsor of this study, put forth the idea that, if the gravity 
potential were known, the same measurements, unsupported by any other information, 
could be used to infer the spacecraft attitude. In May 1993, this idea led to a joint 
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proposal of the Colorado Center for Astrodynamic Research, of the University of Colorado, 
and the Analytical Engineering Co, of Boulder Colorado, to examine the feasibility of the 
suggestion. Subsequently, Goddard Space Flight Center gave a grant to the University of 
Colorado, which in August 1993 issued a subcontract to the Analytical Engineering Co. 
to perform the bulk of the work. This Final Report presents the overall theoretical and 
numerical results of the study. 

The report has 3 main sections. In the 1st, Static Attitude Estimation, the gravity gradient 
tensor is determined in instrument coordinates, as a function of instrument attitude. It is 
shown that, if rotation effects could be somehow removed, spacecraft roll and pitch could 
be determined at the microradian level; but that yaw isn’t observable at all. 

The next main section, Dynamic Attitude Estimation, expands the inquiry to include 
all the rotational effects, but also introduces dynamic estimation, based on the Euler 
equations of rigid body motion, plus the kinematic equations of a spacecraft nominally in 
an observational orbit configuration. The state variables comprise the Euler angles and 
the angular velocity, linearized about this orientation. The Euler equations specifically 
include pitch bias momentum, gravity gradient torque, and torque from unbalanced air 
drag. Air density variations lead to . a process noise, for which a crude power spectrum is 
derived. 

The gradiometer is taken to be composed of an ensemble of accelerometers. For each 
accelerometer, a model is constructed, showing how the output depends on its location, 
the state variables, and the drag. A measurement noise model is given, relating the power 
spectrum to rms noise and averaging time. 

The filter, or estimator based on the plant, noise, and measurement models looks super- 
ficially like a Kalman filter. However, in a dramatic departure from existing theory and 
practice, the filter feedback gains are obtained by minimizing a performance index that 
penalizes the error covariance and the filter settling time. These in turn are computed from 
the power spectra of the various noises. Bryson weighting is used to adjust the penalties 
according to the user’s engineering requirements and desires. This new filter theory has 
previously appeared in print only in [1], which is a condensation of the present report, as 
it then existed. 

Detailed results are presented in Section 3.11. Here it’s sufficient to say that dynamic 
estimation does permit the entire state to be observed, including yaw, for the 3 different 
gradiometer configurations studied; and even in roll and pitch, dramatic gains over static 
estimation are evident. For instance, a 4 single axis accelerometer instrument, with a 
sensitivity of 2 x 10~ s m/s^, averaged over 1 s, could deliver 263, 10, and 29 microradian 
accuracy; in yaw, roll, and pitch respectively; with a filter settling time of 36 s. Worsening 
the accelerometers to 2 x 10 -5 m/s~ still yielded 14, 1, and 2 milliradians; with a filter 
settling time of 63 s. Some 30 cases are examined and interpreted in Section 3.11. 
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The last main section, Self Gravity, examines the disturbing effect of the gravitational 
field arising from the spacecraft’s own mass. Spacecraft components that are fixed in in- 
strument coordinates cause only bias; so only parts that can move are a potential concern. 
Exploratory calculations are made of several common sources of self gravity, including 
articulated devices such as scan platforms, antennas, and rotating solar arrays; thermal 
distortions of the spacecraft structure; and liquids free to move in tanks. 

Articulated devices are shown to yield variable gradient field distortions in the tens to 
hundreds of microradians range; however, modeling, supported by the articulation sensors 
can remove most of the effect. Thermal distortion is typically smaller, but not necessarily 
negligible. Here too, modeling, supported by strain gauges and/or temperature sensors, 
should be adequate to deal with the problem. If liquids are free to move in tanks near the 
gradiometer, they are potentially the worst problem. If nothing is done, field distortions 
of a few hundred microradians are typical. However, with a few extra accelerometers, 
a simple filter can separate the liquid and external fields. Results from a least squares 
analysis are given for 3 different accelerometer configurations, and several liquid— tank 
arrangements. 

The author would like to acknowledge the support of the University of Colorado, partic- 
ularly the Colorado Center for Astrodynamic Research, and its Director, Prof. George 
H. Born, also the Principal Investigator on the Grant to the university on this study. 
Other CCAR personnel contributed substantially, notably graduate student Thomas G. 
Gardner, co-author of [1]; Profs. Penina Axelrad and Donald L. Mackison; and Research 
Associate Dr. Michael E. Parke. Also contributing useful discussions and insights were 
Profs. Daniel B. DeBra and Arthur E. Bryson of Stanford University, Profs. Jason L. 
Speyer and Dino Mingori of U.C.L.A, and Dr. Darrell Zimbleman of Ithaco Inc. And 
there is no forgetting the sponsor, Mr. Seymor Kant of Goddard, originator of the main 
idea, whose continued support throughout, and many hours of consultation, direction, and 
insights were invaluable. 


1.2 Notation 

Uppercase bold roman and greek letters are 2 dimensional arrays, e.g., F. 

Lowercase bold roman or greek letters are column vectors; e.g., r. 

Magnitudes of vectors are non-bold; e.g., r = |r|. 

Lowercase greek subscripts are indices. The Einstein summation convention is used for 
repeated lower case greek indices. 

Overdots signify time derivatives; e.g., x = dxjdt . 
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A T superscript denotes transpose. 

Primes denote scaled variables; e.g., x f . 

Sines and cosines are denoted by s and c respectively. 

A = direction cosine matrix; also, constant matrix in Rjccati equation 

a*, = f c /m = external non-gravitational acceleration on spacecraft 

a* = inertial acceleration of the ith accelerometer 

B = process noise state distribution matrix 

C{ = concern for xt 

C t = settling time concern value 

D(u>) = matrix satisfying Lyapunov equation (78) 

E{x) = expectation of x 

ej — unit vector along axis a in coordinate system x 
F = plant matrix 

f e = external non-gravitational force on spacecraft 

G = matrix defined in (86) 

g(u) — controls vector appearing in (53) 

G = universal gravitational constant = 6.67259 x 10 11 N— m /kg" 
g = gravity field vector at r 

gj — gravitational acceleration at the ith accelerometer 

H — measurement partials matrix 

h = spacecraft pitch momentum bias 

h s — atmospheric scale height 

I n = identity tensor of order n 

J = overall spacecraft inertia tensor 

K = filter feedback gain matrix 

ki -4 — constants defined in (24) 

ks = Boltzmann’s constant = 1.38 x 10 J/K 

kf = air drag force constant defined in (35) 

L — white measurement noise matrix defined in (102) 

1 = set of spacecraft dimensions; also direction cosine vector of disturbing mass 
M = combined “equivalent 75 white noise matrix defined in (98); also a priori covariance 
m — spacecraft mass; also field source mass 
N(a;) = matrix defined in (83) 

P ? — covariance of the error of the estimate 
Q( u) = function of power spectra defined in (74) 
q = filter performance index; also dynamic pressure 
R = measurement error covariance 
R(r) = autocorrelation matrix with delay r 
R(0) = average power in stationary process 
r = field position vector relative to m 
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Ppp = spacecraft center of pressure, relative to center of mass 
r e = earth mean radius = 6.367 x 10 6 m 
S(u>) = general noise power spectrum 
S v , S w = white noise spectra 

s, t superscripts signify spacecraft and trajectory coordinates 
T = absolute temperature 
t = time in seconds 
t' = t/Ct = scaled time 
t s = filter settling time 

U = process noise measurement distribution matrix 
u = vector of controls 

V = bs w v t = white process noise effect matrix 
v 0 = satellite orbital speed 

W = B — KU = process noise effect matrix 
w = process noise vector 

Wd — dimensionless air drag random process’* 

X = F — VM _1 H — linear term matrix in Riccati equation 
x = state vector 
x = estimate of x 
i — dx/dt f 

Y = measurement noise distribution matrix 
y — disturbing mass position vector 

Z = F — KH = observer system matrix 
z — vector of measurements 
6 = error in estimate of liquid location 
€ = variation in spacecraft u; 
r = gravity gradient tensor 

To = Gm/r 3 = gradient scalar due to mass m at distance r 
A = diagonal matrix of eigenvalues A 

A(Z) = eigenvalue of Z; also A = atmospheric density correlation length 
fj, e = Gm e = gravitational constant of the earth = 3.98603 x 10 14 m 3 /s“ 
^ — x — x = error in the state estimate 
a — 3Pi(A) = real part of eigenvalue 
r e = non-gravitational external torque 
T C = measurement noise matrix defined in (88) 

$ = noise matrix integral defined in (90) 

$ = gravitational potential 

ft = process noise matrix defined in (88) 

lo = spacecraft angular velocity 

u) = angular frequency used in power spectra 

u; c = break frequency in power spectrum 

uj h = half power frequency in power spectrum 
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oj 0 = orbital mean motion 


1.3 A Note on Units 

Unless otherwise stated, the units used throughout this report are the SI (International 
System), also known as rationalized MKS units. However, I have also followed common 
practice in the field of gradiometry on the units of the gravity gradient. The natural SI 
unit of gravity gradient is (m/s 2 )/m, or just s~ 2 . Since gradient components at the earth’s 
surface are on the order of 1.5 x 10~ 6 s~ 2 , and are routinely measured to 10“ 9 s~ 2 , or 
better, this has proved to be a rather unwieldy unit. Thus, there has now been world wide 
acceptance of the Eotvos unit 1 : 1 E = 10 9 s 2 . In this report, the SI unit will be used 
everywhere in the formulas; but Eotvos units will be occasionally employed in the text. 


2 Static Attitude Estimation 

i 

2.1 The Tilted Gradient 

Gravitational fields may all be described by a scalar potential field field $. The potential 
due to a particle of mass m at a distance r is: 

$ — —Gm/r (1) 

Note that this potential is negative, but increases toward zero with increasing r. The 
vector gravitational field at this point, due to m, is the acceleration of a free test particle 
there: 

g = -Y$ = -Gmr" 3 r = -T 0 r (2) 

Finally, the gravity gradient tensor field due to m, is: 

r = Vg = To — 13^ ( 3 ) 

The symbol rr^ is an outer product, or tensor product, or dyadic, whatever you re com- 

** rp fy 

fortable with; in contrast to the scalar product r r = r~. 

Outside the earth, the fields are closely approximated by these formulas. Accepting this, 
if the test mass is a spacecraft, in circular orbit about the earth at radius r, then the 

1 Honoring Baron Roland von Eotvos, for his extraordinary experimental work on the equivalence prin- 
ciple in the last century. 
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(4) 


orbital angular velocity, or mean motion u 0 is given by: 

^ = To = Me/r 3 

in which fi e is the gravitational constant of the earth. The actual potential of the earth is 
quite complicated; but differs from (1) by only about 1 part in 1000 in low earth orbit, and 
by even less at higher altitudes. The variations in turn are known to better than 1 part 
in 1000. Thus, if spacecraft attitude is actually inferred from gradiometer measurements, 
this error in knowledge of the field would lead to corresponding attitude determination 
errors on the order of 10 -6 rad, almost surely not the worst error contribution. In any 
case, the intent of the study is to find the accuracy with which a gradiometer can measure 
attitude, given that the field is known. Thus, the study will neglect the effect of field 
knowledge errors. 

On the other hand, neglect of the known deviation from sphericity (mainly oblateness) 
would lead to attitude errors on the order of 10“ 3 radians, usually unacceptable. However, 
the intent of the study is to determine feasibility; so the form of the necessary oblateness 
correction is outside the scope. Later, if feasibility is demonstrated, an add on study, to 
find the form and practical implementation of the correction, would be called for. In the 
same vein, the design of a real system would have to deal with eccentric orbits; but as the 
orbit is not being solved for, the observability of the attitude can’t be seriously affected 
by eccentricity; and the spacecraft orbit will be taken here as circular. 

At this point it’s necessary to introduce coordinates. In general, coordinate systems will 
be described by a set of right handed orthonormal base vectors e£, where a = 1, 2, or 3 
denotes the axis, and x indicates the system. Most important perhaps is the spacecraft 
system e 5 . This is the physical system in the spacecraft to which all the accelerometer 
input axes, and all other instruments, are aligned. For simplicity, it will be assumed that 
the origin of e s is at the spacecraft center of mass. The term spacecraft attitude will 
be taken here to mean the rotation that connects e 5 to a trajectory system e *. In the 
latter system, e* is defined as the local upward vertical, through the origin of e 5 , and 
is parallel to the orbital angular momentum, e£> completes a right handed system, 
and is along the spacecraft velocity vector. It must be emphasized that e* isn’t inertial, 
but rotates uniformly at a rate u> 0 about relative to a system that will be regarded as 
inertial, but won’t need to be identified further. 

The connection between systems may be described by a matrix of direction cosines A: 

= A a 0ep ( 5 ) 

In this study, the spacecraft is assumed to be earth pointing; so A will be taken as a small 
rotation. It then can be expressed in terms of small yaw (^), roll ( 4 > ), an d pitch (0) angles; 
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( 6 ) 


about e^, e!i. and e 3 , respectively. In these terms, and to 1st order in the angles. 

1 6 - 4 > 

A = —0 1 rp 

4 > — 0 1 

From here to' Section 4, the fields will all be due to the earth. The need for e 4 is that g 
and r are most conveniently expressed there: 

g 4 = -IVe^ = IV[-1, 0, Of ( 7 ) 

r t = r 0 diag[2, -1, -1] (8) 

and expressing these in e s , where the instruments reside. 

g s = Ag 4 = ror[-l, 0, -<p] T ( 9 ) 

2 -36> Z<p ' 

r s = Ar 4 A T = r 0 -3* -1 o (io) 

3 <p o -1 

again to 1st order in the angles. These are the transformation rules for contravariant 
vectors and tensors, respectively. Note first, that while pitch and roll turn up in these 
expressions, yaw does not. Physically, this is because r is an axis of symmetry of the fields. 


2.2 Error Analysis 

If we could measure either the gravity or gradient field in e s , we could infer both 6 
and 6. Unfortunately, accelerometers don’t measure gravitational acceleration at all, and 
gradiometers are strongly perturbed by angular velocities and accelerations (see Section 3 
for details). Still, it’s helpful to see how well these angles could be determined if the 
dynamic effects could be removed. For example, if a spot measurement of Tf 3 were possible, 
the error in <p would be: 

6T 0 Sr -i \ 

6<p = — - + 3 <j > — (H) 

ol o r 

Suppose a spacecraft altitude of 500 km. Then r = 6.867 x 10 6 m, and F 0 = 1231 E. A 
gradient component measurement accuracy of .01 E would then contribute 2.708 x 10 
rad to 6<f>. The analysis of 60 leads to the same result, given a measurement of Ffo. In 
each case, the 2nd contribution to the error comes from the uncertainty in the knowledge 
of r. Supposing 6r = 10 m, and <?!> = 0.1 rad, say, this contribution to 6<t> comes to 
4.37 x 10 -7 rad. Since satellite tracking usually leads to determining r much better 
than the horizontal components of position; and attitude control is typically much better 
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than this; the tracking contribution may be regarded as conservative. It will not be 
considered further in this report. Thus, the conclusion of the static analysis is that, if 
spot measurements of the gradient can be made at the .01 E level, then roll and pitch 
determination at the microradian level would be possible. 

If this gradient measurement came from a pair of accelerometers, with an 0.5 m separation 
and independent errors, then their required accuracy would be 

6a = 0.5(10” n )/2 1/2 = 3.536 x 10' 12 m/s 2 

within the capability of the best room temperature accelerometers today, operating in 
space. A pretty stiff requirement; but it will be shown that dynamic estimation allows a 
considerable relaxation. 


3 Dynamic Attitude Estimation 

3.1 Overview 

If gradiometers actually measured the gradient, then a model would be something like 
z = T plus noise, or a subset of its components. A least squares analysis would then 
yield the covariance of the errors in the estimate of the desired <t> and 0, for each discrete 
sample z. However, once it's recognized that any real gradiometer measurement z contains 
functions of and u?, it becomes clear that least squares analysis won’t suffice; and we 
have to resort to dynamic estimation. In the present case, the plant equations take the 
form of the Euler equations of more or less rigid body motion, plus kinematic equations 
relating u> to the attitude angular rates. This structure is developed in Section 3.3 below. 
Actually, as there is very little process noise (external torque variations), these equations 
add considerable strength to the estimates; thus turning a practical necessity into a virtue. 
In the following subsections, these plant equations are developed and linearized, a process 
noise model is spelled out, a filter is synthesized, and the terminal covariance of the errors 
in the estimates is computed; for specific spacecraft, orbit, and instrument combinations. 

A major variation from the earlier gradiometer dynamic estimation studies, [6] and [5], 
is that, instead of treating gradiometers as measuring the intrinsic tensor plus noise, this 
study follows [9] in treating the instrument as an array of accelerometers. The measure- 
ment models then consist of what each accelerometer should measure, plus noise. One 
advantage of this structure is that the measurement noises may now be regarded as un- 
correlated, avoiding the careful treatment needed in [6]. But the big gain comes from the 
much simpler treatment of self gravity, and its detection, to be found in Section 4. This 
model is constructed in Section 3, followed by a measurement noise model. 
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3.2 Spacecraft &; Orbit Models 


For simplicity, the spacecraft will be supposed to be a rectangular parallelopiped, with 
edges l Q aligned along the spacecraft axes e^. Supposing a uniform density p, the space- 
craft mass is: 

m = pl\l<2h (12) 

and the principal moments of inertia are readily shown to be: 

= m(l1 + i3)/12 ; Jo = m{l\ + l1)/\2 ; J 3 = m(/j + io)/12 (13) 

A typical density might be p = 1000 kg/m 3 ; but if articulated solar panels or large antennas 
are present, a lower value would ensue. 

The spacecraft orbit will be assumed circular, at a radius r. For the numerical examples 
an altitude of 500 km will be assumed, for which r =■ 6.867 x 10^ m, ui 0 = .0011095 rad/s, 
and To = 1231 E. No assumption will be made on the orbit inclination, as it doesn’t 
appear in the present analysis. Also, the spacecraft speed in orbit is v c = ru> 0 7614 m/s. 

I 

3.3 Plant Equations 

In [5] it’s shown that the Euler equations of rigid body motion, when modified to include 
an arbitrary bias momentum hw, can be written as: 

Jib = (ju> + hu') x U> + T gg + T e (14) 

in which the external torque has been separated into the gravity gradient torque T gg and 
the nongravitational torque r e , the latter mostly due to air drag. Note that the derivative 
on the left side is the rate of change as seen in e s . Control torques could be included in 
r e ; but as they would then reappear in the filter structure equations, they cancel out in 
the covariance study. 

Unfortunately, this system is nonlinear in u>. Since we are analyzing a nominally earth 
pointing satellite, the nominal value of U) is However, because of the body derivatives, 

a much simpler procedure is to define the variation e by: 

uj — *+- € ( 1 ®) 

Another simplification comes by arguing that, in an earth pointing satellite, bias momen- 
tum, if any, is usually confined to the pitch axis: 

h\y = he 3 ( 16 ) 
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The possibility of additional wheels for control is not precluded by this specification: it s 
only required that their nominal momentum is zero. Substituting these relations into (14), 
and deleting quadratic terms in e, results in: 

Je = o; 0 (Je) x ef + ^(Je^) x (u^e^ + c) + he^ x e + r gg + r e (17) 

Before proceeding with this, it’s helpful to work out t gg . The well known formula in e 
may be put in the form: 

Tgg = 3r 0 e^ x (J'ej) (18) 

Since only J s is readily available, and as what we really need is r s gg , we need to work out 

T s gg = 3r 0 A[ei x (A r J s Ae!)] = 3r 0 (Jn ~ ~ Ju ( 19 ) 

(Jll — J22)0 + J22,4> + >^12 

Note that, while nothing depends on xp, these is a yaw torque, arising from off diagonal 
components of J. These also produce bias torques in roll and pitch. This is why, for earth 
pointing satellites, it is generally preferable to point some principal axis down. Moreover, 
by making this axis (ef) have the least moment of inertia, the gravity gradient torques 
are restoring. In this report, where the main issue is observability, it will be assumed that 
this condition is met. Thus, our assumption for further analysis is: 

J s = diag [ Ji, Ji, J 3 ] ( 20 ) 

In the numerical examples, it will be further assumed that J\ < Jo < J 3 , known to be 
the best configuration for gravity gradient stabilized satellites. With the principal axis 
assumption, the torque reduces to: 

0 

Tgg = 3^ (Jl - (2!) 

^ (ji - j 2 )e _ 

On inserting these expressions into (17), the component Euler equations become. 

Jl^i = [u; 0 (J o — J 3 ) — h]e 2 + r e \ 

J 2 £ 2 = [u; 0 (J 3 — Ji) + h]^l + 3Fo(*7l *“ ^ 3 )^ + r e2 ( 22 ) 

^3^3 = 3To(Jl — J2)^ + r e3 

We can put these in standard form by dividing by the moments of inertia: 

= A:iC 2 + Jf ^1 

k 2 = k 2 C\ + k%4> + J 2 ( 2 ^) 

C3 = k^6 + 
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in which the constants are defined as: 


^1 = 

1 

1 

"o 

JL 

k2 = 

[^£>(^3 — Jl) + h]/J2 

^3 = 

3r 0 (yi - J 2 )/J2 

A: 4 — 

3r 0 (Ji - 


(24) 


J. N U lC Lildt, Wil/Ll LUC UlUOl J i. - ' 

this same vein, if the gravity gradient and other external torque terms are neglected, then 
ei and e 2 decouple from e 3 in (23), resulting in a harmonic oscillator with frequency 

e iven b y ; o (0^\ 

ui'p; = —k\k<2. (25) 

This is the natural nutation frequency, arising mainly from the momentum bias h. 

To complete the plant equations we must add the kinematical relations. With the same 
linearizing assumptions, these are easily shown to be (see for instance [5]): 


ifi = 6l + U) 0 4> 

(p — tn ~ U 0 Tp 

& = €3 


(26) 


We now have a linear system of plant equations of 6th order in e, ip, <t>, and 6. 


3.4 Process Noise Model 

The random process appearing in the Euler equations (23) is the external non-gravitational 
torque r e . At our nominal altitude, this is largely due to air drag; and the random 
component is largely from variations in air density p a . At this writing, actual data on 
high frequency lateral density variations is extremely sparse. The question of a suitable 
air density model for gradiometer studies was addressed in [5]. A flat earth barometric 
model was adopted there: 

p a (r + <5r) = Pa(r)e~ Sr / hs 

where h s is the density scale height. At 500 km, [10] lists p a = 1.905 x 10 kg/m , 
h s — 83 km, and a mean free path of 25 km. These numbers are, admittedly, quite shaky. 
In any case, the dynamic pressure then comes from the speed: 

q = Pa^o / 2 

and with the above numbers, q = 1.106 x 10“ 4 N/m 2 . 
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Since the speed is along eo, and the spacecraft attitude is not far from nominal, the steady 
force from air drag is very nearly; 

fe = ~qhhCD e 2 

Because the mean free path is much larger than the spacecraft, the drag is essentially 
Newtonian, and the drag coefficient is essentially Cd = 2. However, since some inelastic, 
oblique, and diffuse scattering of air molecules is likely, this Cd may be a bit high. In this 
study, the value C D = 1.5 will be adopted, subject to change if a more definite spacecraft 
model is later considered. 

We should also look at radiation pressure. The quantity corresponding to q is / 5 /c, where 
J s = 1360 w/m 2 is the mean insolance at the earth’s mean distance from the sun, and 
c = 2.99776 x 10 s m/s is the speed of light. With these numbers, the “radiation dynamic 
pressure” is 4.54 x 10 -6 N/m 2 . As this value is well below g, and as the variation frequency 
band is much below that for air drag, radiation pressure will be ignored in what follows. 
This point should be revisited if altitudes substantially above 500 km are ever considered. 

[5] goes on to develop a statistical model. It supposes that p a is actually the mean of a 
distribution, to which a random component is added: 

Pr = Pandit) 

where wd(t) is a dimensionless, zero mean, random function of position and time. At 
satellite speed, the spatial variation is much more important. Suppose that Wd(t) has a 
standard deviation a w . What's needed now is a power spectrum. 

Physically, we are looking at dynamic variations in density, with scale lengths on the 
order of h 8 , plus the more or less orbital frequency variation due to solar heating of the 
atmosphere. The latter, while reaching substantial amplitudes, at least for orbits which 
pass close to local noon, is confined to such low frequencies as to have little effect on the 
attitude estimates; so we will ignore it in what follows. As for the dynamic variations, 
we can imagine variability on all length scales, but petering out below distances on the 
order of h 8 - This is the sort of situation that led to the development of the cubic power 
spectrum in [3]; 

s{uj) = ( ! ■ £;) ( 1 + £) (°<"<2^c) (3i) 

and zero otherwise. Here, /£(0) is the average power in the process, and u) c is the break 
frequency where S(uj c ) — S{ 0) / 2. Suppose the autocorrelation of variations falls by half 
at a distance ah s . The time to travel this distance is 

A = ah s /v G (32) 
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( 33 ) 


and [3] shows that, for the cubic spectrum, we should choose; 

7 r 7TV 0 

Wcw = 2A = 2 ah s 

We must also pick 7^(0) and a. For numbers, the best information presently available to 
the author is an analysis of CACTUS data in [11]. Accelerometer data over approximately 
800 s intervals was analyzed at altitudes between 270 and 320 km. Density variations 
of 4%, peak to peak were typical; rising sometimes to ~ 15%, during severe magnetic 
disturbances. The corresponding a w values are .014 and .05. A reasonable balance between 
these values would be ~ .02; but, allowing for a bit greater variability at higher altitudes, 
we will accept a w = .025 as the baseline value. Then, as these time series meet the 
oversampling conditions discussed in Appendix B, /^(O) = — 6.25 x 10 . 

As for a, [11] doesn’t show a power spectrum, but it does show representative time series 
of a normal and a disturbed interval; and it is stated that the apparent wavelengths 
concentrate in the range of 700 to 1500 km. Examination of the time series suggests that 
R(r) falls to 0.5 at r - 50 s. Translating this to 500 km, the corresponding distance is 381 
km, when a - 4.6. Since for a sinusoid, R(r) falls by half at 1/3 of a wavelength, these 
numbers are at least consistent. Again, to allow for a bit more variability at 500 km, we 
will take a = 4 as the baseline. For numerical studies, this leads to = 0.03606 rad/s. 

It remains to convert this to torque. The overall drag force is very nearly; 

f e = —kfl 1 + u>d(^)] e 2 (34) 

where 

kf = pavihhCol^ ( 35 ) 

Supposing a center of pressure at a location r^p in the spacecraft, the torque due this is. 

T e = TcpX f e = fc/(l + Wd(<)]( r cp3, 0, -rcpi] T (36) 

Note that there a deterministic bias force and torque, which must be treated correctly in 
the filter below. Also, while the box structure used here gives rise to no torque along e|, 
an actual spacecraft would likely possess irregularities that would lead to a small torque 
on this axis (propeller torque). To allow for this bit of realism, a component Tqpo can 
replace the zero in (36) above. This provision is included in the filter structure below. 


3.5 Gradiometer Model 

In earlier studies ([6] and [4]), the instrument was modeled as measuring elements of the 
“intrinsic” tensor: 

T — r T u; 2 I 3 — + su) (37) 
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where e is the 3-index permutation symbol. The quadratic oj terms amount to centrifugal 
effects. Because the accelerometers are fixed in e s , there is no coriolis. Here, the instru- 
ment is dissolved into its component accelerometers, partly to avoid the noise correlations 
required in [6] and [4], but mainly to prepare for later studies requiring more elaborate 
error models. This model is also needed in the self gravity studies in Section 4. 

For the purpose of this study, the gradiometer will be taken to be an array of 3 axis 
accelerometers, with their input axes aligned along the e* . For 1 axis accelerometers, it s 
only necessary to select out the appropriate row. It’s convenient to identify a “center” 
of the instrument r c , relative to the origin of e 5 . Then, the ith accelerometer will have a 
position r a i, relative to the center. Thus, its location relative to the center of mass is: 

Ti = v c + r ai (3 8 ) 


Identifying the center in this way is a convenience, useful for entering symmetrical arrays; 
r c should turn out to have no effect on the observability of the state. 

Next, we identify the inertial acceleration. For a perfectly circular orbit, the center of 
mass is subject to -u ;|re t 1 . As for rotation effects, e s is rotating at a rate u>, relative 
to an inertial or non-rotating frame e n . So, purely due to rotation, the inertial velocity 
of the ith accelerometer is (the superscripts indicate the frame in which the derivative is 
observed): 

d n «** --- (39) 


(40) 


fj = — r { ^-r. + wxr^wxrj 
at at 

the latter because r, is invariant in e s . Going to the next derivative 
d n <f n 

rj = — fj = uxri + wx — rj = u> x n + w x (u> x rj) 
dt dt 

Note that u> is the same, whether viewed from e n or e 5 . Finally, on including the external 
non-gravitational acceleration Jig, the ith accelerometer is subject to: 

aj = —u!^re\ -fw x (u> x Fj) + w x rj + a< (^1 ) 


Again, there are no coriolis terms, as rj is fixed in e s . 

On the other hand, the gravitational acceleration of the ith accelerometer is g f plus the 
correction at t* due to the gradient. From (7) and (4), this comes to: 


g i = — u> o re l + 


(42) 


Now, actual accelerometers measure only non-gravitational acceleration; i.e., the difference 
between inertial and gravitational acceleration. These are identical in free fall, when an 
accelerometer measures zero. Conversely, for an accelerometer on a table on earth, it 
measures the acceleration imposed by the table that keeps the instrument from falling 
through the floor. 
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We’re now ready to construct a model of the instrument. For the ith accelerometer, this 
model is: 

Zi = a*-gi + v f (43) 

where v* is the noise in the 3 measurements made by that accelerometer. On substituting 
the above expressions this becomes: 

Zj = x (u; x r*) + u> x r* — Tr* + a e + v* (44) 

Note that the acceleration of the center of mass has dropped out, as would be expected 
for an accelerometer placed there. 

The next step is to linearize this using (15). On neglecting the quadratic terms, and 
recalling that is the same in e n and e 5 , we get: 

z i = UJ 0 (w 0 e5 + e) x (ef x rj) + x (c x r 4 ) + e x n - Tr f + a*. + v< (45) 

We’ll work this out term by term, in the fqjm of matrices of constants times the state 
variables, plus whatever is left over. Starting on the left: 


e-j x (ef x n ) = - r { = -[m, r i2 , 0] 

e x (ef x n) = (rj • e)e%- r< = 


0 0 -rn 

0 0 — r# 

ri i ri 2 0 



Cl 


60 


- 63 - 


ef x (e x rj) = r^e - ezri = 


riz 0 — rji 

0 ri 3 —ri 2 

0 0 0 



J L 


Cl 


60 


_ C 3 _ 


(46) 

(47) 

(48) 


The e term can’t be expressed directly in the state variables; however, by making use of 
the plant equations (23), there follows: 


€ x r i — 


0 ra — r{2 

—riz 0 rn 

ri2 ~ r il 6 


0 fci 0 0 

fco 0 ^3 0 

0 0 0 fc 4 



i — 
f — i 

» 


J i T e \ 


1 

1 

- Ti X 

J^ T e2 
J 3 T e 3 


k 2 riz 0 fc3 r i3 

0 —kiriz 0 

— k 2 rn kyri 2 — fc3 r il 


The T term comes directly from (10): 

2 -30 3 0 

r ri = r 0 | —30 -1 0 

30 0 -1 


~k 4 r i2 
k 4 r n 
0 


n 

Cl 


J 2 l ri3,T e 2 — 


C2 

+ 

Jz^ r il T eZ ~ J\ lr iZ T e 1 

- 

<t> 

0 


JfViOTei — 


(49) 


" 

Til 


ri 3 

-ra 


’ 0 ' 
e 


2 rn 


ri2 

= 3r 0 

0 

~rn 


+ T 0 

-r& 


. m m 


_ rji 

0 




-ri 3 _ 


(50) 
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and on combining all these, and substituting from the process noise model: 


' (fc 2 + w 0 )ri 3 0 —2u 0 r il 0 (fc3-3r 0 )ri3 (3ro - *4^2 

Zi = 0 (uj 0 — ki)ri3 —2u) 0 r i2 0 0 (£4 + 3Fo)ni 

{uj 0 -k2)rn (fci + Uo) r i2 0 0 -(fc 3 + 3r 0 )rii 0 

— 3rji Jn r iZ r cp2 ^3 r i2 r cpl 

+ T 0 0 +k f -J^rarcpi- Ji'rizrcps-m- 1 (1 + Wd(<)l + v* (51) 

rj 3 J\ 1 r i2 r cpZ ~ *^2 r il r cp2 

This completes the description of the accelerometers. There is 1 such 3 vector for each 
accelerometer. 


3.6 Measurement Noise Model 

All accelerometers of the quality needed for this application consist of a case with an 
internal cavity, a proof mass within this cavity, a sensor for determining the relative 
position of the proof mass in the cavity, and a “rebalance” control system for forcing the 
proof mass to the center of the cavity. The force necessaiy for rebalance is the measure 
of the non-gravitational acceleration of the case. Some instruments use a separate proof 
mass for each axis being measured, while others use a single proof mass for all 3 axes. 
The distinction will not be important here, as separate, independent position sensors are 
assumed on each axis. 

There are 3 main sources of noise in this class of instruments — Brownian motion of the 
proof mass relative to the cavity, thermal gradient variations in the cavity, and sensor noise. 
Brownian motion is the result of the thermal energy fcfiT/2 coming to equilibrium with 
the energy stored in the effective spring of the rebalance system. The effect is worsened 
by high temperature T, and by light proof masses. The thermal gradient effect is caused 
by failure of the instrument’s thermal control system to hold down the spatial variation 
in temperature of the cavity. This causes thermal distortion and asymmetry of the cavity, 
leading to errors in the position sensor. A fixed gradient only causes a bias, presumably 
removed by calibration; but variations, either from external heating variations, or from 
noise in the temperature sensors, can be a source of trouble. Finally, sensor noise can 
arise either from the actual sensor amplifiers, or from the power source (either voltage 
or current) used for rebalance. The latter is because this voltage or current is the actual 
measure of acceleration. 

All of these noise sources depend critically on the instrument design. As this part of 
the study is only for feasibility, no specific instrument will be used; but later, if actual 
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instruments are studied, actual power spectra will be needed. Still, we require a power 
spectrum even for a generic instrument; so, for simplicity, a cubic spectrum, of the form 
(31) will be assumed. The discussion in Appendix B shows how the average power K v (0) 
for this spectrum, and the break frequency are determined from the rms acceleration 
error and the averaging time r of the measurement. 


3.7 Filter Structure 

The 1st step in calculating the terminal covariance in a dynamic estimation problem is to 
determine the structure of the filter. This starts with identifying the set of state variables 
that appear in the plant and measurement equations. From (23) and (26), it s clear that 
we should choose: 

x = (ei, 62) e 3> V’) <t>, 0} T 

Following [8], it’s conventional to consolidate the plant equations in standard linearized 
form: 

x = Fx + g(u) + Bw (53) 

Here, F is the plant matrix, u is a vector of controls, g(u), a possibly nonlinear vector 
function, distributes the controls, w is a vector of independent process noises, and B is 
the process noise state distribution matrix. The matrices are readily identified. From (23) 
and (26), we find: 

■ 0 fci 0 0 0 

*2 0 0 0 A '3 

0 0 0 0 0 

F “ 1 0 0 0 Uo 

0 1 0 -u Q o 

0 0 1 0 0 

As for the control and process noise terms, it’s convenient to separate the deterministic 
process noise bias from the random components, and combine them with the actual con- 
trols, if any, to produce the g(u) used here. Since these terms will eventually cancel out 
in the analysis below, the actual controls have no effect on filter performance, and there 
is no need to spell out g(u). Finally, by identifying w with wa(t) in (36), and including 
propeller torque, we have: 

B = kf^cpz/ J\, rqp2 / > — r cpl/*^3) 0) 0] T (^) 

Turning now to the measurement model, the direct appearance of the process noise in each 
of the accelerometer measurements requires a modification of the usual standard model. 

z = Hx + Yv + Uw zb (56) 


0 
0 
k 4 
0 
0 
0 


(54) 
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Here, H is the measurement partials matrix, developed in Section 3.5 above. From (ol). 
this is: 

’ (*T 2 + Wo)fi 3 0 —2 Worn 0 (*3 ~ 3ro)r i3 (3r 0 - fc 4 )r i2 

Hj = 0 (w 0 -kl) r i3 -2u 0 r i2 0 0 (fc4 + 3r 0 )r i i (57) 

(uj 0 — k2)rn (k\ + ui 0 )ri2 0 0 — (kz + 3Fo)rii 0 

and the complete measurement partials matrix is: 

H - [Hf, h£, • • -] T (58) 

For example, if there are 7 vector accelerometers, H will be a 21 x 6 matrix. Again, for a 
1 axis accelerometer, merely select the row in Hi corresponding to the input axis. 

For the measurement noise, it will be assumed that each measured axis of each accelerom- 
eter has a separate independent sensor noise. Thus, v(t) has one element for each element 
of z, and Y is just an identity matrix. A more elaborate model may be found in [2]; so, 
to accommodate specific hardware in later work, Y will be retained in what follows. The 
spectral properties of v(t) were developed aBove. As for the process noise term, having 
established that w is w^t), U comes immediately from (51): 

1 J2 lr i3 r cp2 + Jz ^ r i2 r cpl 

Uj — kf -J^rnrcpi - J^r&rcpZ — m _1 (59) 

Jl^r^rcpZ ~ J<2 l rnrcp2 

The overall U is a column vector with 3 such elements for each accelerometer (or 1 row 
for each 1 axis accelerometer). The remaining terms in (51) constitute the bias zb- As it 
doesn’t affect the covariance analysis below, it need not be spelled out. 

An observer based on these models starts with an estimate x of the state x. This is caused 
to follow the deterministic parts of the plant equations (53), corrected by feeding back the 
residuals, i.e., the actual measurements z minus the measurement model (56). In this case, 
this filter structure takes the form: 

x = Fx + g(u) + K(z — Hx - zb) (60) 

Note that this structure assumes that the control and bias terms are known, and available 
to the filter. The issue buried here is that the controls, whatever they are, are accurately 
modeled in the filter, and that the biases have been accurately determined by some sort 
of in flight calibration technique. Pursuing these points is outside the scope of the study. 

3.8 Terminal Covariance 

The performance of a dynamic filter is generally examined by determining the statistics 
of the error in the estimate, defined by: 

£ = x — x (61) 
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The evolution of £ comes from subtracting (53) from (60): 

$ = Z£ + KYv(t) - Ww(t) 


( 62 ) 


where the observer system matrix is defined as: 

Z = F-KH ( 63 ) 

and the process noise effect matrix is: 

WeB-KU ( 64 ) 

There’s a lot to learn from (62). First, x, x, and all the control and bias terms have 
disappeared, a consequence of the linearizations. Thus, the quality of the estimate doesn t 
depend on the controls in any way, even if they fail to stabilize the plant. This is known m 
the controls business as “the Separation Theorem” , meaning that the design of the filter 
can be separated from the design of the controls (but, alas, not vice versa). The 2nd major 
fallout from (62) is that filter stability requires Z to be stable; i.e., all its eigenvalues are 
in the left half plane. This is a standard requirement in any negative feedback system. 
In filter theory this is put somewhat differently: if a feedback gain K can be found such 
that Z is stable, then the state x is said to be observable by the measurements z. 3rd, 
the diagonal elements and the eigenvalues of Z have the dimensions of inverse time, and 
the filter settling time is essentially given by the inverse of its least negative eigenvalue. 
This fact will be used in the numerical work below to insure that the “optimal filter has 
a reasonable settling time. Finally 4th, since the noises appearing in (62) are all free of 
bias, then after settling, £(t) will also be free of bias. 

Various measures have been proposed to study the quality of the estimate. Here, and 
generally in the references, attention has centered on the covariance of the error: 

P e (t) = £K(t)£ T (t)] ' ( 65 ) 

where E is the expectation operator. The idea that, in a stationary system, Pf(f) would 
have a terminal or asymptotic value, has been around a long time. Calculating this 
terminal value could be quite tedious, if the settling time was long. About 4 years ago, 
William McEneaney, co-author of [6], in unpublished notes, showed that this terminal 
value Pf could be calculated directly from the structural information, and the statistics 
of the noises. After generalizing to allow arbitrary power spectra in the noises, his ideas 
led to [7] and [8], 

The present problem differs from [8] by the inclusion of process noise in the measurement 
model, and various bias terms. Also, [8] required the solution of a set of Lyapunov equa- 
tions involving the autocovariances of all the noises, and it has since been found much 
easier to work with power spectra directly. Since none of these changes and improvements 
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have been documented in any published work, the algorithm for calculating will be 
derived here. 

To begin the analysis, it may be supposed that the filter has been running for all past 
time; so the initial conditions have settled out. Then (62) is solved for “now” in this form: 

£(0) = / e ZM [KYv(^) - Bw(/x)]c?m (66) 

Jo 

Here, the dummy variable p may be interpreted as past time. Strictly, the noise terms 
should be and u; (-/*); but, as we are only after the statistical properties of £, it 

will make no difference. An apparently graver problem is the exponential term — the 
dimensions of tZ depend on those of x, thus calling into question the validity of the formal 
expansion of e tZ . However, from (62), the dimensions of the vector tZx are just those of x. 
Thus, all terms of the form t l Z l x have the same dimensions, and if the exponential is merely 
viewed as a shorthand for the formal expansion, there are no dimensional difficulties. 

The terminal covariance may now be found by substituting this into (65): 

p € = y°° y°° e ZM {KY£[v(M)v r H]Y T K r + W£[w(/x)w T (i/)]W T }e Z ^dfxdu (67) 

This supposes that the expectation and integration operators may be commuted, and 
uses the assumption that w and v are independent and free of bias. On recognizing the 
autocorrelations of the noises, this is: 

P? = l o°° f eZM l KYR ^) YTKT + w] M^) wr ] eZ7Vd / id '' ( 68 ) 


where 


T) = fl-V 


(69) 


Well, autocorrelations and power spectra are Fourier transforms of each other. Using the 
one sided spectra of [7], these relations for any noise component are. 


from which the 


1 [°° 

R ( 77 ) = - / S(u>)c(rju>)du ; 

7T JO 

rco 

S(u>) = 2 / R(T])c(uj7])dr] 

Jo 


average power is: 


1 r°° 

R{ 0) = - / S(u>)duj 
7T JO 


(70) 

(71) 


(72) 


After substituting (70) into (68), and interchanging the order of integrations, there results: 


P £ = - f f f e ZM Q(w)e ZTl/ c(j7u>)<fyr<Wa> 

? TT JO JO JO 


(73) 
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in which: 


Q(u>) = KYS„( W )Y r K T + WS w {u;)W T (74) 

Considerable progress can now be made by a change of coordinates. Combining (69) with 

6 = fi + v (75) 

the double integration region is now the quadrant surrounding the +0 axis, when 

P f = r\[° e Zr > /2 re Ze / 2 Q( W )e ZT ^e- ZT ^ 2 cM)dr ? 

27t Jo — oo J ~v 

+ f e Z,? / 2 [ e Z9/,2 Q (u>)e ZT0 ^d6e~ 2 (inr))dr) div (76) 

JO Jy 

Now, it’s not hard to establish that 

J e Z9/2 Q{uj)e ZTd/2 d6 = 2e Z£l/2 D(w)e ZT< ’ /2 + constant (77) 

where D(a>) satisfies the Lyapunov equation: 

ZD( w ) + DHZ t =Q(w) (78) 

Putting this into (76), and evaluating at the required limits, a considerable simplification 
results: 

p f = -I f°° [° D(w)e -ZT7J c(a.’??)d7j + f e Zv ~D(u>)c(u;T])dri du (79) 

' 7T Jo U-OC' 

and setting rj — » — r] in the 1st integral: 

P f = _I [°° D(w) f e ZTr] c(un i)dr] + [ e Zr] c{ujT))dr 1 'D(uj) du (80) 

^ K Jo [ *70 JO 

when another analytic integral has surfaced: 


leading finally to: 


where: 


[ e’^ jT] c(Lu'T])drf = — (Z + u; 2 Z *) 1 

Jo 


i roc 

P f = - / [NH + N r (o;)ldu; 

7T Jo 


N(w) = (Z + a; 2 Z _1 ) _1 D(o>) (83) 

It may be noted that this analysis would break down in several places but for Z being 
stable. Once again, especially in (81), the dimensions may look flaky. However, letting ui 
represent the dimensions of x iy it is readily shown from the differential equations that the 
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expressions Zyt, Zij/u, and uZ^ 1 all have the dimensions xh/uj. By extension, the ij th 
element of (81) has the dimensions tui/uj. 

This work has established one possible forward procedure. For a given K, Z and W are 
computed from (63) and (64). A set of u values is chosen to cover the region where any 
of the noise spectra are nonzero, with reasonable density. Q(u>) is then determined over 
this set from- (74). Each Q(u/) yields a corresponding D(w) by solution of the Lyapunov 
equation (78), and a corresponding N(w) from (83). Once this is complete, is found 
by integrating (82). 

Since this study began, a remarkable resource has been found to reduce the labor. From 
(82), we can construct the following: 

ZP, + P e Z T = - /°°( ZN(«) + ZN t H + N(w)Z t + N t (w)Z T ) dw (84) 

7T JO 

Next, eliminate N with (83); and, after recognizing that Z and (Z + a;"Z l ) 1 commute, 
the Lyapunov equation (78) reappears. Thus, 

ZP* + PfZ T = G + G r (85) 


G.ifW Z-^QMdu, (86) 

7T JO 

The new procedure calls for computing G directly by (86) from Q(w) and Z, followed by 
solving (85) for Pc. The advantage is that only one Lyapunov equation needs to be solved, 
rather than the 50 or so required for an accurate numerical integration of (82). Note that 
the tempting simplification ZP = G of (85) holds only if Z is symmetric, not the case 
here. 

The main issue now is the integration of (86). This too may be simplified by recognizing 
that Q(u>) is a linear combination of its component noises. Since the power spectra S v (u;) 
and S w (w) are diagonal matrices, with non— zero elements S v k( w) and S w i(u>) respectively, 
we may expand (74) as: 

QH = ^2 SvkW^k + ^ S w i(u>)Qi (87) 

k l 


where 

T k = (KY)ic(KY)[ ; = W/Wf 

in each case meaning the outer product of the indicated column 

G = £* fc T* + 

k i 


( 88 ) 

with itself. With this: 

(89) 
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in which: 


( 90 ) 


<£ = - / (Z+w 2 Z _1 ) _1 5(a;)da; 

7T JO 

The virtue of this dissection is that the integrals need only be taken over the range where 
S(u) is finite. Moreover, in several theoretical cases, analytic solutions for $ are available, 
considerably reducing the arithmetic. 


3.9 Optimal Terminal Covariance 


Having found how to compute from K, we still need to find the feedback gain matrix 
K that yields optimal filter performance. 1st, it’s necessary to quantify “optimal perfor- 
mance”. While P f certainly contains the necessary information, in this 6th order problem 
there are 21 independent matrix elements; so some sort of scalar measure of Pf is needed. 
The software used here is based on a performance index q, constructed from the weighted 
trace of Pf : ■»- 

q = PfrJCl (91) 


In this technique, known as “Bryson weighting”, each C* is the “level of concern” for the 
error For example, if x l were a position, the level of concern might be C* = 1 m. On 
the other hand, Cj = 10 m would show less concern, and cause the optimization to put 
less weight on the variance of Note that the Bryson technique has the virtue that q is 
the sum of dimensionless terms — it doesn’t add apples and oranges. 


A further concern can be added to the performance index q — filter settling time. If 
the K that minimizes (91) also leads to a Z with some eigenvalues whose real parts are 
small (though negative), then we may see from (62) that the settling time of the filter 
will be long, perhaps excessively so. To avoid such a problem, a term may be added to 
(91) penalizing the filter settling time. To see how to do this, consider the behavior of 
the filter evolution equation (62). If A Q symbolizes the eigenvalues of Z, and a a = &(A Q ), 
then the filter response to initial conditions or perturbations may be regarded as a set of 
n exponentially decaying modes, with individual settling times -l/a a . Since all n modes 
decay simultaneously, the overall settling time is given by: 


t s = max(— a,, 1 ) = - (maxcr a ) 


-1 


(92) 


Once again, for the notion of settling to have any meaning, we must rely on Z being stable. 

If overall filter settling is the main issue, we could introduce a concern level C t in seconds 
for the settling time t s , and add t s /C t to q. However, this has the practical effect of causing 
2 or more filter modes to coalesce during minimization. In effect, settling of naturally rapid 
modes is sacrificed to minimize the slowest. I believe a better course is to penalize all the 
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modal settling times equally. To do this, the overall performance index may be taken as: 

< 93 > 

a 

Note: In choosing C u keep in mind that it is the sum of the modal settling times that 
is penalized. In the software, (93) is modified to penalize any <r a > 0 very heavily. This 
adds considerable robustness, but has no effect when Z < 0. 


The added term serves another function. The stability boundary for Z is that all a a < 0. 
Thus, as some a a —< • 0 from the left, t s — » oo. So, adding the C t term erects a barrier in 
the minimizing process against Z going unstable, when > 0 may no longer hold. 

It may be added in passing that an earlier q formulation based on Z -1 instead of Z failed 
badly. The problem was eventually traced to the behavior of complex eigenvalues of Z . 
If A q is real, then as <r Q - 0, the corresponding real part of the eigenvalue of Z 1 (call it 
Pa )-> -oo, as would be expected. However, if A Q is complex, it’s not hard to show that. 
as aa _► 0, p a — ► 0 as well, and doesn’t serve as a proper measure of t s . 

Now suppose we have picked the concerns C, and C u and have a K, such that Z is stable. 
Then P f may be found as detailed above, and q computed. Next, each element of K 
is varied, and the procedure is repeated to get a 6q. Taken together, these constitute a 
gradient of q, relative to the elements of K. A parabolic fit technique is then used to find 
the minimum q along the backward gradient. This whole process is iterated until q has 
reached a stable minimum. The final K is the optimal feedback gain, and the final P* 
represents the filter performance at that gain. 

There is one serious loose end in this procedure — the starting K must yield a stable Z. 
The method initially used in the software is based on Kalman theory. Since S(w) > 0 for 
all u, we can introduce the idea of the half pbwer frequency u h . This is the frequency 
such that half of the total average power is in the interval 0 < w < w/,. From (72), u h is 

given by: , , r*K 

ifl(0) = - £ S(cj)du (94) 

The idea has no meaning for white noise, as f?(0) = oo. For colored noise, S(u>) = 
2R(0)u c /{u>~ + w“), and uj h = u c . For flat bounded noise, i.e., 5(w) = S for 0 < u < 0 
and zero otherwise, u h = fi/2. For cubic noise, from the spectrum (31), it can be shown 
that oj H = 7 u c , where 7 = 0.5327705. Finally, for an arbitrary spectrum, u h is obtained 
from (94), in which ft(0) comes from (72). 

In each case, we form a white noise “equivalent” to S(u). For colored noise, we take the 
level to be S = 5(0) = 2R(0)/u; c . For any other S(u), we substitute the “equivalent 
colored noise — the same fl(0) and w c = u h . Thus, the “equivalent” white noise is flat at 


the level 

5 = 2R(0)/u>h 


(95) 
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Replacing all the noise components with these white noise “equivalents” causes the Lya- 
punov equation (78) to reduce to: 

ZD + DZ r = KYS v Y t K t + WS„W T (96) 

which holds for all u. Since D is now independent of w, P<j may be integrated analytically, 
leading to = -D, when (96) gives a clean connection between K and P ? . On reorga- 
nizing this with the help of (63) and (64), so as to make the dependence on K explicit, we 
have: 

KH P ( + P e H T K r - FPj - P*F r = KMK r - KV T - VK T + BS W B T (97) 


where: 


M = YS v Y r + VS w V q 


(98) 


and: 


V = BS,„U 


(99) 


Since an optimum P* is necessarily stationary relative to variations in K, (97) may be ex- 
pressed in components, and differentiated relative to each K fll/ , leading to this stationarity 
condition for P^: 

KM=P^H r + V (100) 

While this can’t be used directly to eliminate either K or P$ from (97), we need only 

assume that some noise contaminates every measurement component to insure that M is 
non-singular. Thus: . 

K=(P c H r + V)M- 1 ' (101) 

which, except for the V term, is a staple of Kalman theory. When this is substituted back 
into (97), an equally well known algebraic Ricoati equation emerges: 


A + XP C + P^X r = P ? H r M l HP { = P£LP 4 


(102) 


where: 


A s B{S W - S«U r M- l US«)B r 


(103) 


X = F — VM -1 H (1° 4 ) 

All this reduces to the usual Kalman theory when the measurements don’t depend on 
w(t); i.e., U = V = 0. In the software, (102) was solved for P$, and K was- computed 
from (101). While this K is far from optimal for the real power spectra, it does guarantee a 
stable Z to start the real iteration. Indeed, it’s known that, while the Riccati equation has 
many solutions, no more than 1 possesses this property. All this is very interesting; but 
when the performance index (93) was made more robust, it was found that the minimizer 
could find Z < 0 territory, even from a badly indefinite initial guess. 
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3.10 Scaling 


This is quite a large optimization problem. We are minimizing q relative to K. or 
example, if the gradiometer is composed of 4 vector accelerometers, K has 72 elements, 
all of which must be determined. Such problems are touchy, and the difficulties are 
aggravated by poor conditioning in or Z. Some sort of scaling is usually appl ie to 
alleviate this. In the present problem, a natural scaling already exists — the Bryson 
concern levels introduced above. On the hypothesis that the variance P^u is on the same 
order of magnitude as Cf, consider scaling the state variables. 


which are non-dimensional. A similar non-dimensional time may be introduced in the 


same way: 

t' = t/C t 


(106) 


Recall that, in the convention adopted in this report, summation is only over lower case 
greek indices. The covariance of the scaled variables is then: 

Pfo = Eizft) = Ptij/iCiCj) (107) 


The real virtue of such a scaling is that the eigenvalues of should be much closer 
together than those of P£, with a corresponding improvement in the condition number. 
To carry out this scaling, (105) is substituted into (53), leading to: 

*' = F'x' + g'(u) + B'yf (1° 8 ) 


in which ^ /inru 

F'j = FijC t Cjl C{ ; B\j = ; g { = 9iC t /Ci (109) 

It’s not hard to show that the scaling makes all these arrays dimensionless. While it’s not 
necessary to scale the measurements, in the model we must set 


TTx = H' 


moi 


from which 

Hlj = CjHij 


The filter structure then becomes: 


x = F'x' + g'(u) + K'(z - H'x' - z B ) 


in which the derivatives are with respect to t 1 , and 


K'ij = KijC t /Ci 


(HI) 

( 112 ) 


(113) 
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and the error in the estimate 


(114) 


= xj - = it/Ci 

evolves as 

i = Z'i' + K'Yv(t) - W'w(t) (115) 

where 

W' = B'-K'U ; Z'=F'-K'H' (116) 

In components, these matrices are related to the unsealed versions by: 

W' xj = WijC t /Ci ; z' tj = ZijCtCj/Ci (117) 

Note that the matrices Y and U, and thus M aren’t affected by scaling. 

From the determinant relation for eigenvalues, it’s not hard to show that those of Z obey 



A a — CtX a 

(118) 

leading to 


(119) 

and 

t f s = t s j Ct 

(120) 

On substituting these scaling relations into (93), q becomes rather simple. 



q = Tr(P^) + t' s \ 

(121) 


The modified iteration starts by forming B' and F'. Then, transforming the algebraic 
Riccati equation with (107) leads to 

A' + X'P^ + P^X' 7 ’ = P^L'P^ (122) 


in which the primed replacements for A, X, and L are computed as above, except that F, 
B, and H are replaced by their primed equivalents. Note that V — ► V', but no scaling of 
M is required. Solving this leads to a starting value P£ for the main iteration. A similar 
scaling of the u would seem plausible; but a careful review of the formulas shows that it’s 
not necessary. 

Applying the scaling everywhere, the iteration becomes: 

Q'( w ) = K'YS l ,(u,)Y r K' r + W'S^HW' 7 ’ (123) 

The Lyapunov equation is then: 

Z'DV) + D , (w)Z' T = Q V) (124) 
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whose solution leads to N' and P£. In the alternative solution technique, the matrices 
Yfc, Of, <£, and G are merely computed from K', W', and Z'. 

Finally, when q has settled to a minimum, yielding the terminal K' and P£, the unsealed 
covariance is obtained from 

Pzij = CiCjP'a (125) 


and if the gain matrix is needed: 


Kij = CiK'^Ct 


(126) 


3.11 Results & Interpretation 

The calculation of the terminal covariance for a given set of input data requires the exercise 
of 5 programs in sequence, all of which are more or less interactive. These are. (1) 
GRANNY, which asks the user for all the data called for in the above theory, forms the 
various structural matrices, and puts them in the form needed by the later programs; (2) 
ENTRY, which accepts the data package from GRANNY, calculates all the power spectra, 
and constructs the scaled input matrices and arrays needed by the next 2 programs; (3) 
ALRIC, which solves the algebraic Riccati equation (102) to get a starting value for P^; 
(4) TERM, which iteratively computes the feedback gain matrix K which minimizes the 
performance index q, using the algorithm developed in the last 2 sections, resulting in a 
final value of P c ; and (5) GRANPRT, that prints out a summary of all the input data, 
and the results from TERM. Of these programs, GRANNY and GRANPRT are written 
specifically for this study; while the other 3, plus a host of called subprograms, are general 
terminal covariance software. 

During the study, the performance index q, frqm (93), was modified to heavily penalize 
unstable eigenvalues of Z. The new algorithm is contained in a subroutine CQ. This so im- 
proved the robustness of the main algorithm that it was found that TERM could recover 
even from an indefinite Z. Following this, ALRIC was no longer needed. Structurally, 
TERM calls a general function minimization routine DFP, which has several optional 
minimization techniques, including gradient and quasi Newton algorithms. However, nei- 
ther of these performed very well on this problem; mostly because the eigenvalues often 
tend to coalesce near a minimum; when q isn’t a smooth function of K. The problem was 
aggravated here because absolute yaw information comes only from orbital coupling in the 
kinematics equations (26); so, with a reasonable settling time concern value, the settling 
terms in q dominate the covariance terms. 

To get around this, a random jump technique was added to DFP. In this, a controllably 
limited random increment is added to each element of K. A search along the line from the 
old to the new K is made to find a new minimum q, similar to standard gradient methods. 
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The virtue of the random scheme is that the tedious calculation of the gradient is avoided; 
but more steps are needed because there is nothing to motivate the direction of the jump. 

The programs are all written in APL, and implemented on a 486DX 33 Mhz computer. 
Of the lot, TERM (really the called minimizer DFP) is, far and away, the most time 
consuming. It requires only a few seconds for each line search; but anywhere from 10 
to 10 5 iterations are needed to insure a stable minimum. One run required 3 days to 
complete; but most are much shorter. 

Having found a K yielding a “minimum” <?, there is no guarantee that this minimum 
is global. On the other hand, the system designer is assured that the calculated hi er 
performance is achievable; although a different K might, conceivably, yield improvements. 
Another possibility is that the program might have found something like a saddle point; 
so that, a short distance away in K space, q might turn down again. Once something bke 
convergence has set in, a good tactic is to command sets of a few hundred random steps, 
with successive sets at increasingly larger jump sizes. In several cases, a new down slope 
was discovered, occasionally leading to substantially better performance. 

The development of the theory presented here, and its software implementation has been 
long, and, because of the slowness of convergence in many cases, rather frustrating. n 
both theoretical and practical grounds, most important was the discovery that the direct 
calculation of Q(w) by (74) could be replaced by dissolving it into integrals over the 
individual noise components, according to (87 - 90). Another major improvement was the 
discovery that the direct calculation of P 5 by (78), (82), and (83), requiring the solution 
of a Lyapunov equation at each value of w, was much more work than forming G from 
(86), and then solving the single Lyapunov equation (85)i These theoretical improvements 
were eventually combined in the software; so that when all the integrals (90) are analytic, 
numerical integration can be completely avoided. 

About the time these improvements were implemented, an error was detected in the for- 
mulation of the measurement partials matrix H appearing in (51) and (57). While t e 
error might not have caused any serious effect, the results are certainly suspect; and none 
of the 30+ runs made prior to the fix are presented here. 

With occasional variations, the bulk of the cases examined conform to what is called here 
the Baseline configuration. This term covers the satellite orbit, the spacecraft physical 
properties and bias momentum, and the process noise specification. The Baseline orbit is 
circular at 500 km altitude, relative to the mean earth radius. Since the earth is assumed 
to have a spherically symmetric gravity field, nothing depends on the orbit inclination 
The parameters depending on this altitude are w, = .0011094 rad/s, u 0 = 7618.5 m/s, and 
To = 1-2307 x 10" 6 s~ 2 . 

The Baseline spacecraft is taken to be a rectangular box, of dimensions 2.0, 0.7, and 0.5 m 
along the yaw, roll, and pitch axes respectively; and with a mass m - 140 kg. Assuming 
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the mass is uniformly distributed in the box, the principal moments of inertia work out to 
8.6333, 49.583, and 52.383 kg-m 2 . The pitch momentum bias was taken as a rather stilt 1U 
N-m-s, leading to plant constants fci = -1.1587 /s, k 2 = 0.20266 /s, fc 3 - 3.2577 x 10 

Is 2 and fc 4 = -2 8862 x 10 -6 /s 2 . The spacecraft natural nutation frequency is then 
0.48458 rad/s. The spacecraft drag coefficient is 1.5, and the center of pressure offset is 
0.2, .01, and .05 m in the 3 axes, as discussed in Section 3.4. 

The process noise model developed in Section 3.4 led to the numbers p a = 1.905 x 10' 12 
kg/m 3 , k f = 8.2928 x 10 -5 N, air density variation ratio cr w - .025, and the air cell size 
a = 4 scale heights. The resulting process noise break frequency is Ucw = -019204 rad/s. 

The instrument configuration adopted in the early studies consisted of a set of four 3 axis 
accelerometers, located at the corners of a regular tetrahedron. This was done at the time 
the author entertained serious doubts that the entire state would prove to be observable by 
dynamic estimation. Later, when it was shown that as few as 4 single axis accelerometers 
yielded full observability, this configuration (and even richer ones) was abandoned, btill, 
one set of 5 runs was made on the tetrahedral configuration, following all the theoretical 
updates and fixes referred to above; and a table of results will be given. The tetrahedron 
is assumed to be inscribed in a sphere of radius 0.25 m, which gives an edge length of 
0.40825 m. The accelerometer locations r* are then given by the geometric relations lis 
in Appendix A. The instrument center was taken to be at the spacecraft center of mass 
in all cases, although, presumably, this doesn’t matter. 

While the accelerometer noise level was varied between runs, the averaging time was 
always taken as 1 s. The power spectrum was assumed to be cubic, when the discussion 
in Appendix B yields a break frequency of w c = 62.832 rad/s. Unfortunately a fully 
satisfactory analytic treatment of the integral (90) has yet to be developed for the cubic 
spectrum. The existing theory of this integral is presented in Appendix C, and involves 
an eigenvalue-eigenvector decomposition of Z. the problem is that, even near a multiple 
eigenvalue, the program available in APL runs very quickly for eigenvalues alone; but 
when the eigenvectors are also requested, the running time increases dramatically, and 
the program sometimes fails completely. Further research on this integral is intended; but 
there is no guarantee that a fully satisfactory evaluation technique will be found. 

To avoid tedious numerical integrations at each iteration, and get results in a reasonable 
time, it was decided that each cubic spectrum should be replaced by its “equivalent 
colored noise. This is defined as the colored noise spectrum that has the same average 
power R( 0), and the same half power frequency u> h \ i.e., the frequency within which ha 
of the total power is found. For colored noise, the solution of (94) is - u; c ; while lor 
cubic noise, if we let u h = ~fv c , then the substitution of (31) into (94) leads to 


y 4 — 47 3 + 167 — 8 = 0 


(127) 


whose only positive real solution is 7 = 0.532770504326. Thus, if we have a cubic spectrum 
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for which w c is given, the “equivalent” colored noise has a break frequency given by 7 <*>c- 
Lacking a better approach, an option for doing this has been temporarily included m 
ENTRY. The key advantage of the colored noise substitution is that a complete analytic 
evaluation of (90) is available when Z is definite. A sketch of the proof is given in Appendix 
D. The result is included in the routine CQ for evaluating q. 

Results for the tetrahedral configuration are given in Table 1 below. In all cases, the 
concern level for the angular velocity u was taken as 10' 6 rad/s, tighter I now think than 
most applications require. As in most places in this report, the concern level for roll and 
pitch was 5 x 10' 5 rad/s. On the other hand, recognizing that yaw is both harder to 
measure and less critical in many applications, the yaw concern was taken as .005 rad/s. 
Except for Case 5, the settling time concern of 10 s was adopted. The reader should eep 
in mind that this concern is applied to the sum of the settling times of all 6 modes of the 
filter; and, typically, 3 or 4 dominant modes would have similar settling times. 

In all the tables, the 1st column is the Case number, sometimes with a (), indicating a 
note below. The 2nd column is the standard deviation of the measurement noise. Column 
3 contains the final value of the performance index q. Column 4 provides the settling time 
of the slowest mode. Columns 5-7 contain the standard deviations of the components of 
co, while Columns 8-10 give the corresponding results for the angles. 

Now for some discussion of the results. The accelerometer for Case 1 represents a very 
high quality instrument, if not the absolute best (and most expensive). Ordinary inertial 
grade accelerometers, if space qualified, probably couldn’t deliver this performance, even 
with the greatly lowered scale factor errors from operating in free falL Nevertheless, on 
examining the performance, the settling time is found to be quite unsatisfactory. Moreover, 
on dividing the 6 following standard deviations by their corresponding concern levels and 
summing, it may be seen that q depends on settling time, almost to the exclusion o 
anything else. The message seems clear — low’settling times are intrinsically difficult to 
achieve. This conclusion is certainly in line with the observation that i> is only observable 
through roll-yaw coupling, whose characteristic time is l/w 0 = 901 s. However, read on. 

Case 2 raises the measurement noise by an order of magnitude. Sure enough, the error 
levels, with the exception of a q, all increase; but they still fail to contribute significally 
to g; and both t 3 and q increase moderately. Case 3 raises measurement noise by another 
order of magnitude, but entering the range where accelerometer costs might not dominate 
the attitude control system. To my intense surprise, both q and t s dropped, although, 
except for u> 2 , all the errors worsened; and q continues to be dominated by settling times. 
This case caused a great deal of soul searching. After a lot of checking failed to find any 
further errors in either the theory or the software, it was concluded that the performance 
seen in this table is achievable; but it probably doesn’t represent the best; i.e., global 

minima for q. 
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Continuing to Case 4, measurement noise was again raised by an order of magnitude, 
representing the crudest instruments likely to see this kind of service (and the lowest cos 
i/this casef the run was interrupted, printed out with GRANNY and later resum«L It s 
instructive to list both results in the table, and they are referred to as Cases 4A an • 
This time, everything but i> has worsened, and the settling time dominance contmu . 
On continuation, a further reduction of q by 3.6% was observed to a result that was 
stable enough to insure at least a local minimum. Various small adjustments are > seen in 
the errors, while t, was dropping by only 1.9%. An examination of the eigenvalues of Z 
showed that settling is dominated by 2 close complex pairs, and the improvement between 
4A and 4B was only 2% in the worst pair; but 5.7% in the slightly better pair. I believe 
that the message here is that an overall improvement in performance can sometimes e 
achieved, at the cost of worsening some error components that aren t contributing 

to g. 

The final run, Case 5, returned to the noise level of Case 3; and, in an attempt to improve 
settling, lowered C t by an order of magnitude to 1 s. The rather start ng resu wa 
success of a sort, in that t. dropped by nearly half, at the cost of having all the errors nse 
dramatically, except for o; 3 which had been unusually high in Case 3. Settling still domi- 
nates q, but not by so overwhelming a margin. Because of the change in C t , the mcreas 
in q over Case 3 shouldn’t be viewed as a worsening of performance. An examination of 
the eigenvalues shows that settling is dominated by 2 complex pairs m Case , 
shifted to a single complex pair plus a single real in Case 5. This was the 1st appearan 
of what might be a different valley in K space. 

On reviewing these findings with the sponsor, he felt that this exploration should be 
abandoned, in favor of studies of sparser, and lower cost instruments. At this tun , 
the features for making the q calculation more robust had not yet been added so 1 g 
random jumps in K space were hazardous, in %at they might lead to an indefinite Z, and 
what amounted to a program crash. Thus, although clearly desirable, a search for better 
tetrahedral gradiometer solutions was never carried out. 


Table 1 — Tetrahedral gradiometer 


Case # 
O=note 

BBS 

X 

<? 

<3 U1C X 

t s 

s 

nrad/s 

- — o 

u>2 

nrad/s 

nrad/s 

rj) 

/xrad 

4> 

^xrad 

6 

H rad 

i 

0.2 

244.1 

786.2 

13.24 

0.5723 

9.479 

1.836 

0.8203 

0.7304 

0 

2 

371.3 

928.9 

5.851 

2.522 

50.74 

61.6 

27.91 

12.97 

3 

20 

230.0 

624.1 

8.912 

0.2826 

122.3 

216.1 

48.78 ; 

20.79 

O 

4A 

200 

464.9 

1173 

30.1 

3.274 

128.4 

104.1 

69.81 

32.89 

4B 

200 

448.1 

1151 

30.58 

3.376 

137.3 

93.0 

66.15 

35. 15 

5 (1) 

20 

1081 

351.1 

72.75 

6.934 

2.281 

2933 

243.1 

57.15 


(1) Ct lowered to 1 s 
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In deciding on candidate sparse configurations of accelerometers, we were guided by (10), 
which tells us that the 12, 13, 31, and 32 components carry useful information; but that the 
latter 2 only repeat the former. Of course, this is because T is symmetric. This suggests 
that the minimum configuration could take either of 2 forms. In the 1st, called here the 
“yaw gradiometer”, a pair of roll input accelerometers, separated in yaw, is augmented by 
a pair of pitch input accelerometers, also separated in yaw. The 2nd configuration, called 
here the “cross gradiometer”, consists of 4 yaw input accelerometers, 2 separated m roll, 
and the other pair in pitch. There are other possible arrangements of 4 accelerometers, 
but they haven’t been explored. 

In this static view, with similar separations and sensitivities, these 2 configurations should 
yield the same performance. However, as shown in the next section, they don’t act the 
same in the face of self gravity disturbances. Moreover, in dynamic estimation, since the 
intrinsic tensor (37) isn’t symmetric, equal performance of the 2 configurations shouldn’t 
be expected. 

In going from Table 1 to Table 2, the principal change is the switch to the yaw gradiometer. 
The Baseline orbit, spacecraft, and process noise assumptions are all the same, as are the 
angular and angular velocity concern levels. However, except for Case 9, all the settling 
time concerns are set to 1 second. For comparison, the only entry in Table 1 with this 
settling time concern is Case 5. The overall result was nothing short of startling — the 
drop in settling times by more than an order of magnitude is very hard to understand. 

Some insight comes from from looking at the physical arrangement. In the yaw gradiome- 
ter, the separation was taken as 0.5 m, along the yaw axis. From the tables in Appendix 
A, the corresponding separation in the tetrahedral configuration is only 0.28868 m. On 
the other hand, the tetrahedral configuration has 12 measurements, vs. 4 in the yaw gra- 
diometer. However, from the standpoint of the gradient tensor, only 8 of these make any 
direct contribution to observing the angles; so the “overkill” is only 2 to 1, relative to the 
yaw gradiometer, for a relative advantage of %/2. On combining these factors, the yaw 
gradiometer appears to have a relative advantage of 1.2247, hardly enough to explain the 
improvement in performance. 

In the author’s view, the main difference is that K has 72 elements in the tetrahedral 
configuration, vs. only 24 for the yaw gradiometer; and we are searching for a minimum of 
q in a space of this many dimensions. Clearly, the opportunities for additional local minima 
will be much greater in the former case. In support of this view, note that the error in u >3 
in Case 5 is improbably low, suggesting that the solution there is far from the best. Also, 
an examination of the eigenvalues (not shown in the tables) indicates that the solutions 
in cases with identical assumptions tend to fall into groups. The interpretation that each 
group represents a separate valley is hard to escape. Unfortunately, the robustness features 
were still to come when the runs of Table 2 were made. A lot of additional insight might 
result from reruns of a couple of cases from Tables 1 and 2. 
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Cases 6-8 were the 1st test of the yaw gradiometer, and differ only in the accelerometer 
noise levels. The surprisingly good performance in Case 8, with the most noise, suggests 
that the results in Cases 6 and 7 are suboptimal; and indeed, the eigenvalue constellations 
look completely different. With the trial assumption that Case 8 is near optimal, let s do 
some analysis. 

1st, the roll and pitch performance are better than that predicted by the static theory m 
Section 2.2 by a factor of about 18,000. Can dynamic estimation really yield so large an 
improvement? Well, yes. Suppose there were no process noise. Then (74) shows that by 
merely choosing K = 0, we would get Q(u») = 0, when P € = 0; i.e., the terminal estimate 
is perfect, in spite of noisy measurements. This is actually a common occurance m filters 
that compute K on the fly. If, through unwarranted optimism, the filter decides that 
is small, it will compute a correspondingly small K. It will then tend to ignore the 
measurements, and “go to sleep”. This is essentially what happens in the current theory, 
if the effect of the measurement noise dominates that of the process noise. 

1st, consider pitch. From (36), the torque noise is a w k frcpl = 4.146 x 10“ 7 N-m. The 
pitch measurement noise comes from the pair of roll input accelerometers; and expressing 
this as torque, the noise is V2J 3 a v /6l = 2.963 x KT 5 N-m, where SI = 0.5 m is the 
separation along yaw. However, the filter also knows that the measurement noise has a 
break frequency of 33.475 rad/s, vs. uj cw = .019204 rad/s for the process noise, whose ratio 
is 1743. The filter can then safely conclude that the bulk of the accelerometer output is 
high frequency noise, rather than the result of process noise, and choose gains to suppress 
it. Assuming it could do this perfectly, its estimate of the residual measurement noise 
would be (Tar = 2.963 x lO" 5 /vTf43 = 7.097 x 10 -7 N-m, comparable with the process 
noise. Then, with the assumed pitch dynamics, the corresponding angular error would be 

erg ~ o'cr/{^Z UJ ^w) — 3-674 x 10 5 rad 


within a factor of 7 of the “optimal” solution. Not bad for an arm waving argument. 

The argument for roll is a bit more subtle. This time, the dynamics come from (23). On 
removing the slow gravity gradient term, the roll dynamics are 


o k2 T el . , 

6 2 -f r “ ~ K f 


J 1 


J o 


^2*1 u d (t) + r -f-w d (t) 
J i *>2 


the last from (36), with the constant term removed. Since w d {t) ~ w OT w,i(i), the 2nd term 
is negligible with the Case 8 numbers, and the particular integral of the solution may be 

written as t 

e 2 (t) = f 2rc ? ~ j w d (T)s[uN(t -T)\dr 

Jl^AT Jo 

With a page of algebra, the statistics could be worked out in terms of S m (u;), and evaluated, 
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so: 


but this is only a rough argument, 





kfk‘rr C p?,'u>d{t) 


kfk or^cr,., (8.293 x 1Q- 5 )(0.2027)(.05)(.025) _ 9.628 x 10 6 rad 

a * ~ yfijyuiN^ \/2(8. 633) (0. 4846) (.0192) 2 

The corresponding argument for the roll measurement noise is similar to pitch. The 
angular acceleration noise level is the same as for pitch: y/2a v /6l. However, the dynamics 
are now those of nutation; so the roll standard deviation, as inferred directly from e 

measurements, is 

_ v/ ^° v _ x 10 — l = 2.409 x 10" 6 rad 
* “ u%6l (0.4846)(0.5) 

but once again, the filter will know that the bulk of this is high frequency measurement 
noise, rather than the result of process noise; so the same argument as for pitch yields 
as/V 1743 = 5.77 x 10 -8 rad. Of course, the optimization process wouldn t waste the 
filter resources on driving < 7 * so low; but the argument certainly shows why relatively 
crude accelerometers should be able to deliver microradian performance, and that static 
estimates have little to do with the performance to be expected from dynamic estimation, 
at least when there isn’t much process noise. 

The other question is how the filter can settle in 20, s, when the characteristic time for 
yaw determination is 901 s? To make this more plausible, recall that settling time doesn t 
pertain to acquisition, a nonlinear process, but rather ,the recovery of the filter from a 
small unmodeled perturbation. The size of the perturbation doesn’t matter, so long as 
the linearization assumptions in the theory aren’t violated. So settling is the time for the 
filter to pile up enough information to reduce Its step error by 1/e. This information is 
accreted from the measurements, and depleted by the process noise. 

Here, (26) tells us that measurements of 4> ought to yield information about at orbital 
rate. Thus, if we could ignore the process noise, in 1 settling time we should achieve 
a performance H ~ ^/(o/ 0 t a ); and with the numbers from Case 8 this is .0003 rad, 
compared to .0008 rad in the actual run. Since the process noise for this Case has been 
shown to have a slightly greater effect on the filter than the measurement noise, our arm 
waving argument appears to be quite good, as is the argument that Case 8 is nearly 
optimal. In general, if we are satisfied with crude performance in yaw, compared with 
roll, then we can get in a time that is short compared with the orbital period. On the 
other hand, if we demand equivalent performance, we’ll have to wait for something like an 

orbital period to get it. 

In Case 9, the settling time concern was weakened to 10 s, relative to Case 8. This 
gave the disturbing result that t s slightly improved, while all the errors increased, the 
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opposite of what was intended. The only reasonable conclusion is that the result is seriously 
subop timal, a conclusion backed up by the drastically different eigenstructure. 

Case 10 was the same as Case 6, except that the process noise <r w was raised by a factor of 
4. Physically, this was an increase in the variability of air density to a quite unreasonable 
level. The result was that everything turned sour, including the settling time, which 
increased by a factor of 5. In the light of the settling time discussion above, this makes 
excellent sense, in that information about yaw is being removed about as fast as it s being 
created. This result also suggests that at satellite altitudes well below 500 km, steps should 
be taken to streamline the spacecraft, or improve its symmetry, so as to reduce r^. 

Cases 11-20 started out as a run to settle the question of whether moving the gradio meter 
away from the spacecraft center of mass would make any difference. All these runs followed 
Case 6, with the difference that r c = [1.0, 0.5, 0.4] r m, instead of zero. When Case 11 
differed substantially from Case 6, in spite of the belief that r c should make no difference, 
it was decided to repeat the run. Since these cases didn’t consume much time, 10 such 
cases were run in all, without clarifying the original question at all, but strengthening the 
notion of multiple minima considerably. Note in particular that Cases 15 - 18 do appear 
to lie in the same valley; but then, why were they consecutive? At the request of the 
sponsor, the inquiry was given up to adopt a set of parameters closer to current spacecraft 
needs. 


Case # 
()=note 
6 

7 

8 

9 (1) 

10 ( 2 ) 

11 (3) 

12 (3) 

13 (3) 

14 (3) 

15 (3) 

16 (3) 

17 (3) 

18 (3) 

19 (3) 

20 (3) 


Table 2 — Early yaw gradiometer results 


<T y 

nm/s 2 

q 

ts 

S 

nrad/s 

LJO 

nrad/s 

prad/s 

jrrad 

<t> 

fii&d 

9 

H rad 

20 

74.75 

23.24 

27.93 

22.88 

1104 

1674 

10.6 

4.443 

2 

60.14 

12.01 

49.64 

29.14 

835.1 | 

1261 

12.91 

5.243 

200 

106.4 

19.58 

158.7 

46.49 

315.8 

802.8 

6.587 

5.484 

200 

9.025 

16.64 

335.8 

138.5 

493.9 

1229 

14.28 

11.47 

20 

267.6 

114.2 

458.1 

340.3 

1670 

17,621 

71.65 

26.85 

20 

57.43 

15.67 

36.08 

21.01 

158.0 

454.4 

2.603 

4.899 

20 

31.7 

7.92 

183.4 

24.31 

434.7 

1380 

11.33 

6.201 

20 

51.97 

11.23 

72.79 

4.746 

133.5 

116.9 

5.558 

6.026 

20 

81.5 

19.78 

94.52 

23.03 

115.2 

690.9 

9.725 

4.455 

20 

49.19 

9.812 

79.35 

35.39 

267.4 

610.4 

9.34 

3.505 

20 

54.12 

9.739 

78.53 

35.37 

263.9 

612.9 

9.295 

3.471 

20 

48.89 

9.804 

78.61 

35.37 

264.7 

612.1 

9.298 

3.475 

20 

48.81 

9.751 

78.11 

35.27 

263.3 

610.9 

9.252 

3.477 

20 

36.52 

9.037 

38.02 

25.03 

119.4 

291.1 

3.658 

4.339 

20 

21.8 

5.806 

76.69 

24.04 

477.8 

917.8 

4.682 

5.133 


(1) Settling time concern = 10 s 



(2) Raise process noise a w to 0.1 

(3) Offset instrument to [1, 0.5, 0.4] m from [0, 0, 0] 

In going to tlie final set of runs, several major changes were made. 1st, the robustness 
improvements in CQ were completed, and the use of ALRIC to start the iteration was 
abandoned. 2nd, detailed notes were taken during each run, showing the progress m 
reducing q, and the maximum random step size and number of steps in each automatic 
sequence. This is an interactive feature of the minimizer DFP, in which the user is asked to 
provide this information; the program carries out the commands, and shows the running 
progress on the screen. On completion of a sequence, it pauses to permit a review of the 
current results. During these final runs, a technique was developed of alternately using 
large steps to locate new valleys, and small steps to refine the results. Once q appears 
to have stabilized, several thousand steps should be taken, with increasingly larger steps, 
until the running screen results show frequent jumps into the unstable region (indefinite 
Z). If this shows no further improvement in q, the current minimum is likely to be global. 

The 3rd major change was in the concern levels. After a lengthy consultation with the 
sponsor, it was decided that the earlier runs had far too much concern for angular velocity 
errors. Accordingly, the concern was weakened by 3 orders of magnitude to .001 ra /s 
in all 3 axes for all the remaining runs. However, it was recognized that some mission 
applications would require tighter angular rates, and thus lower concern values. In most 
applications, yaw determination isn’t as important as roll and pitch. The values chosen 
here for most of the final runs was .001 rad in yaw, somewhat tighter than before, and 
5 x 10 -5 rad in roll and pitch, as in most of the earlier runs. However, the angle concerns 
were loosened by an order of magnitude for Cases 26 and 27. The settling time concern 
was taken as 10 s for all of the remaining runs. 

The 4th and last major change was that the line search routine, called by DFP during 
each step, was completely rebuilt. Now, if a jump failed to improve, the opposite direction 
is also tried. Then, if either improves, successively larger jumps are made along the line, 
until q turns upward. A parabolic fit is then made to the last 3 points, and q is again 
computed. Whichever actual q calculation yields the least value is kept as the result of 
the line search. Large jumps rarely yield improvements; but when they do, the drop can 
be dramatic. 

For all the final runs, the Baseline orbit and spacecraft were used. The atmospheric 
model was also the same, except for increased process noise in Case 28. The pitch bias 
momentum was Baseline, except for lower values in Cases 29 and 30. The instrument was 
a yaw gradiometer, except for Case 23, in which the cross configuration was examined 
Various accelerometer sensitivities were employed; but the averaging time was 1 s in a 
cases. The lowest measurement noise level, 2 x 10 -8 m/s 2 , is not the best available, ut 
was regarded as the most expensive instrument likely to be used in attitude determination. 
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Cases 21 and 22 were actually the same run, featuring the best accelerometer. In Case 21, 
about 2000 steps of varying sizes brought q from the unstable region to about q = 52, and 
another 2000 to 49. Then, another 1000 at a large step dropped it to 28. Another 1000 
medium and small steps yielded a stable q = 2G.2, the listed result in Table 3. Because 
of a mistake, the program had to be restarted as Case 22, leading almost immediately to 
a new valley, with q = 24.5; and after another 2000 medium and small steps, stabilized 
at q — 18. 1G. Unfortunately, no more large steps were attempted, to insure that this 
minimum is global. Moreover, there is nothing quite like this case to compare to, as the 
similar runs in Table 2 were all made with C* = 1 s. Except for settling time, the results 
are quite impressive, since all the errors are well below their concern levels. Even a settling 
time of 36 s isn’t too bad; and it suggests that a bettef solution might have been found. 

Cases 24 - 27 successively raised the measurement noise, to see how fast the performance 
degraded as the instrument was cheapened. At the highest poise levels, Cases 26 and 
27, lowered expectations caused us to weaken all the angular concern levels by an order 
of magnitude. Except for an initial drop in t Sy further impugning Case 22, it steadily 
degraded, until Case 27, where it jumped to an impractically high level. By and large, 
except for a minor glitch at Case 25, and a major one at Case 27, the errors worsened 
monotonically. While it’s unlikely that the Case 27 result is global, an improvement in t s 
would almost surely be accompanied by a further worsening of the errors. Thus, Case 26, 
with a measurement error of 2 x 10 “ 5 m/s 2 , and yielding milliradian angular error levels, 
is probably the minimum instrument quality that should be considered. 

Case 28 differed from Case 22 by raising the process noise by a factor of 4, similar to Case 
10. Nothing dramatic happened to the errors, and t s even dropped, lending weight to the 
argument that Case 22 wasn’t a global minimum, although no large jumps were tried here 
at the end either. Nevertheless, it’s likely that Case 28 did find the global minimum. 

A quite different experiment was tried in Case 29. The pitch momentum bias was lowered 
from 10 N-m-s to only .08 N-m-s, relative to Case 22, in order to reduce the nutation 
frequency ujjq to .005 rad/s, only 4.5 times u> 0 . The result was disastrous, in that t s rose 
to 3887 s, although the errors were certainly tolerable. Since some 3500 large steps were 
made at the end of this run, without effect, the result is fairly credible. What appears to 
have happened here is that, in reducing the stiffening effect of the momentum bias, the 
roll - yaw behavior becomes much more susceptible to process noise. The message is that, 
if we wish to get away with so low a bias, either a much higher altitude, or a much cleaner 
spacecraft, or both would be required. 

Because of the huge change from Case 22 to Case 29, a final Case 30 was thrown in with 
the intermediate value of 1 N— m— s for the bias. This proved to be an extremely tedious 
run, involving at least 60,000 medium steps, and several days. And it didn t appear to be 
quite done then. This sluggishness discouraged any attempt to find another valley with 
large steps, especially as the cause of this behavior isn’t known. As for the results, they 
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are intermediate between Cases 22 and 29, lending much credibility to the conclusions 
reached from Case 29. 


Table 3 — Final runs 


Case # 
()=note 

nm/s 2 

9 

^5 

S 

nrad/s 

o/2 

nrad/s 

a>3 

nrad/ s 

rp 

fi rad 

4> 

/xrad 

e 

jzrad 

21 

20 

26.21 

48.55 

759.9 

176.1 

732.5 

60.19 

24.29 

70.76 

22 

20 

18.16 

36.33 

170.4 

550.1 

1561 

262.5 

10.11 

28.87 

23 (1) 

20 

29.47 

43.98 

5777 

20,118 

50,059 

1369 

184.3 

42.9 

24 

200 

15.15 

24.08 

1479 

2612 

10,966 

1928 

29.45 

73.04 

25 

1000 

30.15 

40.23 

17,447 

2695 

3997 

323.0 

50.49 

110.7 

26 (2) 

20,000 

48.54 

62.64 

366,788 

58,463 

86,995 

13,607 

1028 

1828 

27 (2) 

200,000 

555.1 

1253 

238,845 

99,403 

9316 

7117 

2398 

2615 

28 (3) 

20 

11.27 

24.2 

284.0 

1104 

2286" 

224.0 

16.37 

25.05 

29 (4) 

20 

812.0 

3887 

358.6 

1489 

2957 

77.53 

24.3 

6.274 

30 (5) 

20 

52.15 

103.7 

644.2 

539.7 

1302 

335.8 

39.12 

27.83 


(1) Cross configuration 

(2) Angular concern levels raised to [.01, .0005, .0005] rad 

(3) <j w raised to 0.1 

(4) Bias momentum reduced to .08 N-m s 

(5) Bias momentum mid value of 1.0 N-m-s 


4 Self Gravity 

4.1 Discussion 


In Section 3.5, a model of the gradiometer was developed, showing how variations in the 
local gravitational acceleration field would afTect the individual accelerometers. However, 
there it was tacitly assumed that the field is entirely due to the planet. In fact, the 
spacecraft generates its own field, with its own substantial gradient. To get an idea of 
the importance of this, consider a sphere of radius r, and density p. Then, from (2), the 
gravitational scalar for this sphere at its surface is: 


r 0 = 


Gm 

^ 3 " = 


G 47T 3 

-3 T r P = 

r 6 3 



(128) 


and is independent of the radius. Thus, if a bowling ball had the same density as the 
earth, the gradient tensor due to the ball would have the same value at its surface as the 
field at the surface of the earth, or in low earth orbit. (A real bowling ball is less by a 
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factor of about 2.) It should be clear from this that, while typical spacecraft densities run 
1 or 2 orders of magnitude less than the earth, its internal field can seriously distort the 
external field, and lead to unacceptable attitude determination errors. In general, in tfie 
inertial instrument field, this effect is known as “self gravity”. 


The simplest, and most conservative way to look at this problem is to ask how much the 
field is tilted by a spacecraft component, and thus how bad the error would be if the 
component were completely ignored. Suppose the component is a point mass m, at a 
location y relative to the center of the instrument. Then from (3), the disturbing gradient 
is 

/ 3yy r 


6T = r n 


-4 


where 


r 0 = Gm/y 3 


( 129 ) 


( 130 ) 


Next, from the discussion in Section 2.2, the roll error that would arise from ignoring this 
disturbance would be 

(131) 

3 t Oe t Oe 

where Toe is the gravitational scalar of the earth, and 


U = Vitv 


( 132 ) 


Similarly, the corresponding pitch error would he 

(133) 

l 0 « 

Thus, if is located exactly on any instrument axis, there is no error; but if it were 1 
in out, and at equal angles to the 3 axes, then li — 1 / n/ 3 ; and at 500 km altitude, 1 kg 
would cause errors of 181 microrad in roll and pitch. Clearly, just about any spacecraft 
design would contain components that cause unacceptable errors if ignored. Fortunately, 
there are several measures available to the system designer for alleviating the problem. 


By and large, for the spacecraft under consideration here, the spacecraft mass is fixed in 
instrument coordinates, and thus causes a set of accelerometer biases that are indistin- 
guishable from other instrument biases. Overall, these biases need to be determined by 
some sort of in flight calibration system; but the possible design and performance of such 
a system are not part of the present study. 

Of present concern are spacecraft components that are free to move relative to the instru- 
ment; e.g., solar panels, scan platforms, and articulated antennas. Other possibilities are 
propellants or other liquids that are free to move about in tanks, and thermal distortions 
of the structure. In most of these cases, some sort of modeling is possible. For instance, 
the field of a solar panel can be computed from a model of its structure; this together with 
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the output of its shaft angle encoder yields the acceleration disturbance at the location 
of each accelerometer. While no such modeling is perfect, the bulk of each acceleration 
disturbance can be removed in this way. 

4.2 Point M ass Displacement 

To deal with self gravity in some generality, suppose a mass m is at a location y in e s 
(instrument coordinates). Then, the gradient tensor at the origin, due to m, is given by 
(129). Now suppose that m is displaced by a small vector 6y. The change 6T in the 
gradient can be worked out as 

6T = T(y -1- Sy) - T(y) = ~ \ySy T + Syy T -I (l 3 - 4 yyT ) y ^ y ( 134 > 

y L \ v* / 

to 1st order in 6 y. 

The information on the roll angle is, from (10), 

3roe4> = ri3 (135) 

Thus, the apparent change in <f> y due to 6y is: 

6<f> = ^ 1^3 + h e i ~ 5 / 1 / 31 ^ 6 ^ Fo/Toe = /F o/Toc (136) 

in which 1 is again the direction cosine vector corresponding to y, and € is a sort of strain 
vector corresponding to 6 y: 

e = 6y/y (137) 

What are the worst combinations of 1 and €? Well 1 is a unit vector; but for the question 
to make sense, it’s also necessary to limit e. Probably the easiest way to do this is to 
introduce a limiting isotropic strain e > 0: 

€ T € = c 2 (138) 

One way to proceed is to assume that, temporarily, 1 is a given unit vector, and find the 
values of € obeying (138) that lead to extreme values of /. The variational Hamiltonian 
for this formulation is: 

7i = / 1 ^3 + + A(€^c — c 2 ) (139) 

and the necessary conditions for a stationary / are: 

dH/dcx = h - 5/^3 T 2Aci = 0 (140) 

dH/dei = -511*2*3 + 2Ae 2 = 0 (HI) 
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(142) 


dHIdtz = li - 6Iil| + 2Ac 3 = 0 

On using the 2nd condition to eliminate A from the 1st and 3rd: 



1 

CM 

II 

cs 

to 

i) £2 

(143) 


5l2 ; 3«3 = (*><3 - - 

l) £2 

(144) 

and on substituting these into (138): 



WM 2 = 

l§ (5 1? - l) 2 1 (5/i/ 2 i3) 2 I i? (5<1 - l) 2 «2 


— 

25l?l§ (if + If 1 If) 

- 20 $1 -1 l\ + ll\ e\ 


- l 

(r.ifil + 1 ? 1 if) <f 

J 

(145) 

whose solutions are 




e 2 = ± r ohhh( (W?ls + l l + l l) 

, -1/2 

) =±51i hhQ 

(146) 

leading to 

El = ±*3 (51? _ 

1 )<? 

(147) 


e 3 = ±h ( 5 I 3 - 

l)Q 

(148) 


Proceeding from these we can construct 

l T e = ±3 hhQ (149) 

and with a little more algebra we get 


/ = Te 2 /Q = T£\/ 5, l / 3 + Z 1 + l l ( 15 °) 

Evidently there are 2 opposite values of e leading to opposite, but equally bad extremes 
of /. 

We can now ask what values of 1 yield the worst of the worst, given that e is so demonically 
picked? Well, by inspection, it may be seen that the most painful choices are h = 0; 
/j, f 3 = ± l/v/2 independently; for all of which: 

/ = T3e/2 ; € = ±el (151) 


and the worst angular distortions are: 



The geometrical interpretation is now straightforward. If either the yaw-pitch or the 
pitch— yaw components of T are measured to determine <fi y then the worst places to put m 
are in the roll plane, at equal distances from the axes, when the worst strain is radial, i.e.,. 
toward or away from the origin. 

The corresponding analysis of 60 is based on Fi2i and clearly leads to the same worst cases 
as for <f> y except that / 2 and I 3 are interchanged, and the worst problems are when r?z is in 
the pitch plane. We may also conclude that no single 1 is worst for both 0 and 4 > . Overall, 
it appears reasonable to conclude that, for arbitrary locations and displacements, it would 
appear prudent to assume variations of order 

|M, \60\ ~ tTo/Voe ' ( 153 ) 

For example, if the nominal 500 km altitude is assumed, then Toe = 1-231 x 10” s ; 
and if m = 10 kg, and is 1 m from the gradiometer, then To =-6.673 x 10 10 s 2 ; when a 
0.1 m displacement yields field distortions of order 5.42 x 10 rad. A scan platform on 
a small satellite would be in this range; while loose propellant in a tank might be a few 
times worse. 

A simple application of these ideas is thermal expansion of the spacecraft structure. Sup- 
pose the spacecraft is unevenly heated on one side, causing 100 kg of the overall mass 
to be displaced radially from the instrument. If the structure is aluminum, with a co- 
efficient of thermal expansion of 2.5 x 10” 6 /I<; and the temperature variation is 20 K, 
then c = 5 x 10" 4 . If the 100 kg is 0.5 m from the gradiometer, To = 5.34 x 10“ 8 s” , 
and the field distortions are about 2.17 x 10” 5 rad. If this level of error can’t be safely 
ignored, then a thermal model of the structure, supported by thermocouples, or simple 
strain gauges, or both can be used to calculate corrections at each accelerometer. 

Whenever (153) yields a result that is too large to be ignored, the designer may resort to 
modeling. That is, if 8y is measured, the full nonlinear model of g(y) may be computed 
at each accelerometer location, and used as a correction. A reduction of the error by a 
factor of 10 - 100 should then be possible, depending on the accuracy of the model, and 
of the measurement of 6y. 

An issue that has been buried here is that the field distortions might have a different effect 
in dynamic rather than static estimation, partly because the angles reappear in the plant 
model, and because there are additional state variables to be estimated. It’s arguable 
that, since the only source of information on the angles comes from I\ the field distortions 
should be carried through to the estimates of 0 and 0 without much change; although 
this says nothing about 0 or o>. The matter clearly needs testing by further analysis or 

simulation. 
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4.3 Rotating Dipole 


For many spacecraft configurations, the largest articulated mass will be a solar panel. If 
the panel were modeled as a flat plate, rotatable about an axis lying in the plate, then 
the field of the plate could be found by integrating (3) over the surface of the plate. This 
field could then be evaluated at each accelerometer location, as a function of the rotation 
angle. While these accelerations might be substantial, the designer might reasonably 
hope that the variations in the field with rotation would be smaller, especially if the 
rotation axis passes through the center of the plate. As a practical matter, the analytic 
expressions for the field components are pretty complicated; but a table of corrections for 
each accelerometer and rotation angle would only amount to a few thousand numbers. 

To get a handle on whether this modeling is necessary, the problem has been greatly 
simplified. The gravitational efTect of the plate has been modeled by condensing the plate 
mass into a pair of point masses m, separated by a massless roll of length 2 g, orthogonal 
to, and centered on the rotation axis. Clearly, if g = 0 there is no variation. Suppose 
the rod center is at a location r c = [u, v, in) 7 , relative to the center of the gradiometer. 
Let the variable rotation angle be a. Finally, for definiteness, suppose the rotation axis is 
parallel to pitch, i.e., e|. Then, the position vectors of the masses are. 


r + = [u 4 gca, v \ gsa, w] 7 

(154) 

i T 

r_ = [u — < 7 ca, v — r/sa, w| 

(155) 

and the magnitudes of these are 


r+ = + 2g(uca + t>sa) -(- g 2. 

(156) 

r?. = r\ — 2g(iicot + usa) T g 2 

(157) 


The gradient tensor due to this mass dipole is the sum of 2 terms of the form (129). Thus 
the total distortion in 6 is: 


r 12 Cm \ (u + qca)(v + gsa) (u - gca)(v - gsa) 

<5# ~ ~ rr - — — r 175 " ,- 5 

31 Oe 1 Oe [ 1 + - 

and the corresponding total distortion in <p is: 



r i 3 Gmw fu-\-gca _ u - gca\ /ir;Q) 

6<t> - TTTT- = -p I 15 + -5 I V ' 

While 69 and 6<j> possess analytic derivatives with respect to a, even Mathematica couldn’t 
solve for the stationary points analytically. However it was easy to plot both functions vs 
a, and several cases were examined. The extreme values of the distortion were taken from 
each run, and their difference computed. The results are given in Table 4 below. All the 
dimensions are in meters, and the distortions are in microradians. In all cases m = 5 kg, 
corresponding to a 10 kg solar panel. For other values of m, just apply the appropriate 

ratio to the table. 
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Table 4 — Simplified solar panel model 


Case 

u 

V 

W 

g 

^^rnin 

^max 

A0 

^$min 

50max 

A<t> 

i 

1 

1 

1 

0.5 

24.9 

35.4 

10.5 

26.2 

49.2 

23.0 

2 

1 

1 

0 

0.5 

62.7 

196.4 

133.7 

0 

0 

0 

3 

0 

1 

1 

0.5 

-12.4 

12.4 

24.8 

-27.1 

27.1 

54.2 

4 

l 

0 

1 

0.5 

-12.4 

12.4 

24.8 

71.4 

98.9 

27.5 

5 

l 

l 

0.5 

0.5 

48.0 

108.5 

60.5 

51.2 

161.1 

109.9 

fi 

0.5 

0.5 

0.5 

0.3 

171.9 

278.1 

96.2 

368.2 

881.9 

513.7 

7 

0.5 

0.5 

0.5 

0.5 

49.1 

177.2 

128.1 

98 

1262 

1164 

8 

0 

0 

1 

0.5 

-47.6 

47.6 

95.2 

0 

0 

0 


Of course, in a real case, the full integration over the actual panel shape should be carried 
out, and g at each accelerometer, rather than T computed; but the table should give 
a pretty good idea of the size of the distortions to be expected. As a general conclu- 
sion, except in extreme cases such as 6 and 7, distortions on the order of 0.1 mrad are 
to be expected; reducible by an order of magnitude or more by careful modeling. On 
the other hand, the extreme cases tell us that a design that permits moving masses too 
close to the instrument can cause either unacceptable errors, or overly stringent modeling 
requirements. 


4.4 Loose Liquids 

If the spacecraft design includes 1 or more tanks, partly full of liquid propellants or 
cryogens, then migrations of the liquid can dominate the self gravity concerns. Lacking 
any physical constraints other than the tank walls, an accurate model of the motion is 
hard to imagine; and sensors that might locate the liquid aren’t readily available. I his 
problem was faced during a JPL study of a Lunar Orbiter, that proposed to measure 
gravity with a gradiometer; and employed a spare Mars Observer spacecraft, containing 
some rather large propellant tanks. It was immediately recognized that the self gravity of 
the propellant was a very serious concern. 

Since the gradiometer was of the differencing accelerometer type, the author suggested 
that, with a few extra accelerometers, it might be possible to estimate the external gradient 
and the location of the propellant simultaneously. The subsequent analysis led to |9), in 
which it was shown that the separation is possible for a wide range of instrument designs. 
The problem here is actually simpler, in that the external gradient is assumed known, but 
the attitude is not. Here, the basic idea of [9] will be followed; but the problem will need 
to be reformulated. 

The liquid will be modeled as a spherical blob of known mass m, and at a location y 
relative to the center of the instrument. If the liquid doesn’t wet the tank walls, then 
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m > 0. Otherwise, it plates out on the walls, leaving a spherical ullage space in an 
otherwise full tank. The problem is the same, except that m < 0, and the instrument bias 
will include a contribution from a completely full tank. The estimation state would then- 
include the unknown y, and the attitude angles appearing in the accelerometer model. 
Actually, in the interest of linearization, we may suppose we have an estimate y of y, and 
an error in the estimate: 

6 — y — y (160) 

Prior to any measurements, we may take y as the center of the tank. In these terms, we 
may take the estimation state as: 

x = \(t>, 0, <5 1 , <5 2) . <53] T (iG!) 

As this is a static analysis, only the static terms in the accelerometer model (44) will be 
retained, plus the liquid disturbance: * 

Zt = — r e rt 1- W{ 1 Vj (1G2) 

where T e is the external gradient, and w f is the liquid disturbance field at r<. As before, 
Vf is the accelerometer measurement error; but here it’s a sample error, rather than a 
random process. 

The 1st term comes from (10): 

— r e r; — r<te|Mr i2 — 3</>rj3 - 2r*i, 30rfi I r* 2 , H3 ~ (163) 


Next, the disturbance at the ith accelerometer is 


w i = Cm 


|y - Hi 3 


On letting 


hi = y- r { 


this may be written as 


Wj = 6*771 


hi -6 
I hi - 6| 3 


and on expanding to 1st order in 6: 

Wj = Tj [hi + 3(lf«)U-6] 


where 

and 


Tj = Ginh i 3 
li S h i/hi 


(164) 

(165) 

(166) 

(167) 

(168) 
(169) 
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All is now ready to construct the measurement partials matrix H. There is 1 row in H 
for each 1 axis accelerometer; 

H 4 (j) = dzi/dx (17Q), 

where j is the input axis number for the ith accelerometer. From the above relations, this 
works out to: 


Hj(l) = 3 — To e riZ, r 0e r i2? 1 

Ht(2) = 3 0, Foe r tli F ri^? 2 - , I 

H<(3) = 3 -roeni, 0, Tikila, F^a, r <(&“^)] 


(171) 

(172) 

(173) 


Suppose the a priori knowledge of the angles is pretty crude; say cr<^ — ag — 0. 1 rad. Also, 
“somewhere in the tank” corresponds to around a y = .05 m in each axis, for an 0.3 m 
diameter tank. Overall, this gives an a priori covariance of the error in x of 


M=diag[.01, .01, .0025, .0025, .0025] 


After a least squares analysis has been performed on the measurements z, maximum 
likelihood theory shows that the a posteriori covariance of the error in x is given by: 

P- 1 = M' I |H T R' 1 H (174) 


where R is the covariance of the measurement errors v. This is the information form 
of a covariance update, in which the N4 term is the a priori information on the state x, 
and the remaining term is the information contributed by the measurements. It should 
be pointed out that maximum likelihood theory, on which this formula is based, requires 
that the measurement errors v have Gaussian distributions, not generally the case. More- 
over, Bryson weighting of P, as in the dynamic estimation theory, would probably yield 
lower angular errors, especially as the liquid location errors aren t intrinsically interesting. 
Unfortunately, the theory for this has not yet been worked out. 

The theory of (171 - 174) has been incorporated in a computer program LIQ. It interac- 
tively enters the satellite altitude, the location and input axis number of each accelerom- 
eter, the standard deviation of each accelerometer error, the location of the center of the 
tank, and the mass of the free liquid blob. Since the accelerometer errors are independent, 
R is diagonal, and made up of the variances of the measurement errors. The program 
computes H and then P, and lists the angular and location standard deviations. In addi- 
tion, if m were at y, and ignored, the angular errors would be given by (131) and (133). 
Of course, these errors are independent of the measurement error. All these results are 
given in the tables below. 
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Of the many possible instrument configurations, 3 were studied in some detail. The 1st 
of these, called here the full yaw gradiometer, consists of a pair of 3 axis accelerometers, 
mounted at 2 locations along the spacecraft yaw axis. The separation was 0.5 m for alb 
cases studied. The results for this configuration may be found in Table 5. Some discussion 
is needed to make sense of the results. Column 1 gives the case number for each run 
of LIQ, to identify it below in the text. The parenthetical numbers refer to notes below 
Table 4, but used in all the tables. Column 2 contains accelerometer measurement error. 
Column 3 is the location vector of the center of the tank, in meters. Columns 4 and 5 
are the post measurement standard deviations of the roll and pitch errors respectively, in 
inicroradians. Columns C - 8 are also standard deviations, and show how well the liquid 
blob has been located. Finally, Columns 9 and 10 are from (131) and (133) respectively, 
and show the errors that would ensue if m were at' the tank center, and no modeling were 
performed. Except as noted, all runs assumed m = 100 kg, and an altitude of 500 km. 

The 1st series of runs, Cases 101 - 107, is a sequence of increasing measurement noise. 
They show that measurements better than 1(T 7 m/s 2 are needed to extract the angles 
better than the a priori estimates; and better than 10 9 m/s 2 if 1 mrad performance is 
required. This is also the level at which locating m begins to become important, as shown 
both from the consequence of ignoring m entirely (0.348 mrad), and by the location errors 
beginning to drop below their a priori values. 


It should be noted here that a value of 10 -9 m/s 2 in Column 2 doesn’t mean that the 
instruments must perform this well. If, say, the accelerometers deliver independent mea- 
surements every second, then they could be averaged over the time needed for significant 
motion of the blob. After disturbances have settled, this might be 100 sec in orbit, when 
this performance could be delivered by instruments only capable of 10 m/s . Such 
thinking suggests that a better job of removing this type of self gravity error might be 
done by dynamic estimation, based on a fluid mechanics plant model; but the idea won’t 
be pursued here. 


Next, Cases 108 and 109 varied the satellite altitude to 1000 and 300 km respectively, 
compared to Case 104 at 500 km. As might be expected, there is a more or less linear 
variation in the recovered angles, and in the efTect of ignoring m completely. On the 
other hand, there was almost no effect on the recovery of the blob location; although the 
measurement noise was not low enough to give significant improvement over the a priori 
value. 


Cases 110 - 112 reduced m to 10 kg, relative to Cases 102 - 104. The effect was exactly 
an order of magnitude reduct ion in the error from failing to model m at all, just as would 
be expected. At the same time, the location error was seriously worsened, especially for 
better measurements — smaller masses are harder to find. As for the angle recoveries, 
there was only a small (1% ) improvement, and that only for the best measurement. Case 
113 further reduced in to zero, relative to Case 103, with consistent results. 
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Case 114 reduced the tank center distance by half, relative to Case 103, and in the same 
direction. The error of neglect increased by a factor of 8, as should be pretty obvious. 
Finding the blob becomes a great deal easier, as would be expected; but, curiously, the. 
angular errors worsen by only about 10% — although the disturbance is much worse, the 
filter has an easier time removing it. 

Finally, Cases 115 - 120 varied the tank center location relative to Case 104. In general, 
the angular recoveries worsened by a few per cent, mainly because the tank locations 
were all a bit closer to the instrument. Note that, in Cases 115 and 118 - 120, the 
disturbance is ignorable, as predicted by (131) and (133), indicating that there are some 
preferred tank locations. Alternatively, if tanks can’t be put in these desireable locations, 
the gradiometer orientation could be changed, although linearization of the angles would 
then be more complicated. 


Table 5 — Full yaw gradiometer 


Case # 
()=note 

<?V 

nm/s 1 2 3 

y 

m 

6<f> 

/xrad 

60 

/Trad 

6x 

mm 

6y 

mm 

6z 

mm 

/xrad 

101 

.005 

i, i, i 

3.976 

4.007 

2.204 

1.747 

2.148 

347.8 

102 

.01 

i, i, i 

7.951 

8.012 

4.393 

3.488 

4.282 

347.8 

103 

0.1 

i, i, i 

78.66 

78.82 

32.22 

28.6 

31.7 

347.8 

104 

1 

i, i, i 

766.7 

766.7 

49.6 

49.49 

49.58 

347.8 

105 

10 

i, i, i 

7638 

7638 

50 

49.99 

50 

347.8 

106 

100 

i, i, i 

60,814 

60,814 

50 

50 

50 

347.8 

107 

1000 

i, i, i 

99,159 

99,159 

50 

50 

50 

347.8 

108 (1) 

1 

i, i, i 

946.7 

946.7 

49.6 

49.49 

49.58 

429.4 

109 (1) 

1 

i, i, i 

701.7 

701.6 

49.6 

49.49 

49.58 

318.3 

110 (2) 

.01 

i, i, i 

7.866 

7.882 

32.22 

28.6 

31.7 

34.78 

111 (2) 

0.1 

i, i, i 

76.68 

76.67 

49.6 1 

49.49 1 

49.58 

34.78 

112 (2) 

1 

i, i, i 

766.1 

766.1 

50 

49.99 

50 

34.78 

113 (3) 

0.1 

i, i, i 

76.61 

76.61 

50 

50 

50 

0 

114 

0.1 

.5, .5, .5 

87.08 

87.85 

5.228 

4.196 

4.85 

2782 

115 

1 

0, 1, 1 

767.3 

767.3 

49.48 

48.48 

48.51 

0 

116 

1 

1, 0, 1 

772.3 

767.5 

48.34 

49.22 

47.94 

958.4 

117 

1 

1, 1, 0 

767.5 

769.0 

48.31 

48.1 

49.22 

0 

118 

1 

1, 0, 0 

814.8 

814.8 

26.29 

41.35 

41.35 

0 

119 

1 

0, 1, 0 

766.1 

796.9 

47.12 

47.82 

45.92 

0 

120 

1 

0, 0, 1 

796.9 

766.1 

47.12 

45.92 

38.24 

0 


&@ign 
H rad 


347.8 

347.8 

347.8 

347.8 

347.8 

347.8 

347.8 

429.4 

318.3 

34.78 

34.78 

34.78 

0 

2782 

0 

0 

958.4 
0 

0 

0 


(1) Cases 108 and 109 are at 1000 and 300 km altitude, respectively. 

(2) Cases 110 - 112, 208 have m = 10 kg. 

(3) Cases 113, 209, 307 have m = 0. 
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The 2nd configuration studied varies from the 1st by moving the 2 accelerometers with 
yaw input axes to the roll axis; the idea being that a 2 dimensional configuration might do 
a better job of separating local effects than the 1 dimensional full yaw gradiometer. The. 
arrangement is now a cross, with 0.5 m arms. 'I he results are given in Table 6. 

The 1st series, Cases 201 - 207, are a repeat of Cases 101 - 107. The performance is 
minutely better in roll, and significantly better in pitch; although, curiously, the blob 
location has slightly degraded. The reason is fairly clear from (10) the yaw input 
axis accelerometers yield no angular information when separated only in yaw; but when 
their separation is in roll they give pitch information. The possibility of putting these 
accelerometers instead on the pitch axis wasn’t studied; presumably, the benefit would 
switch from pitch to roll. 

Case 208 is a repeat of Case 112, in which m = 10 kg. As in Cases 201 - 207, the only 
change is an improvement in pitch, presumably for the same^reason. Case 209 further 
reduced m to zero, with no change in filter performance. The reason is evident the 
error from neglecting 10 kg is already too small to make much difference. Case 210 repeated 
Case 114, in which the distance y was halved relative to Cases 203 and 103. The results 
followed the same basic trend: the improvement in location accuracy keeps the angular 
recovery from getting much worse. Finally, Cases 211 - 216 repeat Cases 115 - 120, in 
which various locations y were examined. There were no surprises; the location accuracies 
were more or less the same, if a bit scrambled; and the improvement in pitch recovery was 
again demonstrated in all 6 runs. 
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The 3rd, and last, configuration studied is called here the triple cross. It’s an extension of 
the cross configuration for which dynamic estimation results were obtained above. There, 
it consisted of 4 accelerometers with yaw input axes, arranged in a 0.5 m cross on the roll 
and pitch axes. Here, a 3rd yaw axis arm is added, again with yaw input axes, and 0.5 m 
separation, making a triple cross. The results are in Table 7. 

Cases 301 - 306 again repeated Cases 101 - 106. For accelerometer accuracy of 10“ 9 m/s 2 
or worse, there was essentially no difference, since location isn’t improved, and the extra 
accelerometers don’t contribute angular information. For more sensitive accelerometers, 
the location improvement was not as good as the full yaw gradiometer, and this in turn 
worsened the angle recovery. 

Case 307 is like Case 304, except that m = 0. As in Case 209, the reduction in the 
disturbance was too small to have much effect on the angles. Case 308, like Cases 210 
and 114, halved y. Once again, the dramatic improvement in location accuracy kept the 
angles from getting much worse; and again the triple cross proved inferior to either of the 
other configurations. Finally, Cases 309 - 311 mirrored Cases 118 - 120 and Cases 214 - 
216, in which various instrument axis values of y were tried. Ihe results showed a modest 
scatter, with no obvious message beyond the observations above. 
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Table 7 — Triple cross of yaw accelerometers 


Case # 

<?V 


y 


6<f> 

<50 

6x 


bz 


Ovign 

()=note 

nm/s 2 


m 


/irad 

jxrad 

mm 

mm 

mm 

jxrad 

jirad . 

301 

.005 

1, 

i, 

1 

7.491 

7.491 

7.633 

16.56 

16.56 

347.8 

347.8 

302 

.01 

1, 

i, 

1 

12.8 

12.8 

14.76 

25.75 

25.75 

347.8 

347.8 

303 

0.1 

1, 

i, 

1 

79.37 

79.37 

47.57 

38.06 

38.06 

347.8 

347.8 

304 

1 

1, 

i, 

1 

766.7 

766.7 

49.97 

49.4 

49.4 

347.8 

347.8 

305 

10 

1, 

i, 

1 

7638 

7638 

50 

49.99 

49.99 

347.8 

347.8 

306 

100 

1, 

i, 

1 

60,814 

60,814 

50 

50 

50 

347.8 

347.8 

307 (3) 

1 

1, 

i, 

1 

766.1 

766.1 

50 

50 

50 

0 

0 

308 

0.1 

.5, 

.5, 

, .5 

116.4 

116.4 

8.458- 

11.47 

11.47 

2782 

2782 

309 

1 

1, 

o, 

0 

800.7 

800.7 

.'22.7 

50 

50 

0 

0 

310 

1 

o, 

1, 

0 

766. 1 

806.5 

37.54 

47.84 

50 

0 

0 

311 

1 

o, 

0, 

1 

806.5 

766.1 

37.54 

50 „ 

47.84 

0 

0 


A Appendix: Polyhedron Formulas 

Early in this study, the gradiometer was composed of accelerometers located at the vertices 
of a regular polyhedron. Tetrahedra, cubes, and octahedra were to be considered, but 
only the tetrahedron was actually tried. The program GRANNY, which enters all the 
problem description data, then required that the user enter the instrument geometric 
data, including the number of faces n of the polyhedron (4, 6, or 8), and the instrument 
size. The latter is given by either the polyhedron edge length /, or the radius of its 
circumscribed sphere r, or the polyhedron volume v. 

If r or u are entered, l is calculated by a set of formulas taken from Mathematica, and are 
summarized in the following table: 


71 

4 

6 

8 

l/r 

yg/3 

2/\/3 

v/2 

uW] 


1 

V 2/9" 


These formulas are built into GRANNY. 

As for the locations of the vertices, they are assumed centered on the origin according to 
the following tables, where each column is a location vector: 


Tetrahedron: 


l 

J5 


11-1-1 

l-i-i l 
l-i l-i 
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Cube: 


1 

2 


Octahedron: 


11 1 1 -1 -1 -1 -1 
11 - 1-11 1 - 1-1 
1 - 11 - 11 - 11-1 


10 0-10 0 
0 0 1 0 0 -1 
010 0 -1 0 


l 

71 


When a central accelerometer is added, a zero columo is appended to the appropriate 
matrix. If the instrument isn’t at the center of mass of the spacecraft, an offset vector is 
added to each column of (this possibly augmented) matrix. 

B Appendix: Averaged Measurement Noise 


The instruments studied in this report are modeled as measuring the acceleration of their 
case, plus random noise. In practice however, they generally average the analog output for 
some length of time r, and deliver a digital result at the end of each interval. The feasibility 
study considers only analog instruments, and thus takes r = 0. On the other hand, the 
instrument manufacturers often characterize their devices as delivering ‘ samples (really 
averages) every r seconds, or alternatively, at a sample rate of 1/t Hz. The noise associated 
with these averages is then specified by a standard deviation cr. This appendix deals with 
relating this type of specification to the parameters of the assumed cubic power spectrum. 

This situation was examined in [12], where it was found that for an arbitrary noise power 
spectrum S(w), the variance of the averages is given by: 

a 2 = — — rj f S(w)[l - c(rw)]^ (175) 

nr 1 J o w 

On assuming the cubic spectrum (31) for the analog noise the result can be put in the 
form: 

a 2 — R(0)f s { TU! c) (176) 


where: 


«■» !jf(* 


du 

2 


(177) 


This integral may be evaluated analytically in terms of the sine integral function: 
/,(*) = |si( 2 ») I- | (“ I cx) - 1(1 + s 2 x) 


(178) 



in which: 


, , ( v sz , y , y 

Sl ( y ) - Jo ~* dZ V 3 • 3! + 5 • 5! 

The function looks ghastly for x « 1; but it actually behaves quite well: 


( 179 ) 


/.<«)= l-£ + °(« 4 ) (180) 

This is the oversampling limit; i.e., if a time series is very frequently measured, but is long 
enough to cover many cycles of the highest noise frequency, then 71 ( 0) is the variance of 
the samples, and the distinction between sample and average disappears. Actually, this 
limit holds for any S(w), as is readily seen from (175). ^ 


The other limit, x » 1 is also pretty clean: 

Si(.x) -> tt/2 ; f s {x) -* tt/x ^ ( 181 ) 

Overall, f s {x) is a monotonic decreasing function, whose behavior can be seen from the 
table: 


X 

0 

0.1 

0.2 

0.5 

1 

2 

3 

/s( X ) 

1 

0.99956 

0.99823 

0.98901 

0.9574 

0.84917 

0.71822 

X 

5 

10 

20 

50 

100 

200 

500 

/s( x ) 

0.50907 

0.28422 

0.14958 

.061632 

.031116 

.015633 

.006271 


When a was measured by the manufacturer, the repetition frequency 1/r was probably 
chosen about an order of magnitude below the break frequency w c /{2-x). Adapting this 

reasoning, we can pick: , DO , 

20 tt/t ( 182 ) 


LO r 


so that TU, C = 207T = 62.832 rad; and 71(0) = .0492401a 2 . This assumed structure has 
been used to determine the measurement noise power spectrum in the study. 


C Appendix: Cubic Noise Integral 


I’ve not been able to find a direct analytic solution for the matrix integral $ in (90). 
However, suppose Z is subject to an eigensystem decomposition: 

Z = EAE -1 ( 183 ) 

where A is a diagonal matrix of the eigenvalues A j of Z, and E is a matrix whose columns 

are the corresponding eigenvectors. This is always possible unless there are repeated eigen- 
values, a quite unlikely case. However, complex eigenvalues and vectors have frequently 

occurred in the study. Now, x 

Z- 1 = EA _1 E -1 ( 184 ) 
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and we have: 

3> = — E f (A f - u; 2 A~ l )~ l S(u)dojlS~ 1 

* Jo 

Since the integrand is clearly diagonal, this becomes: 

3> — Ediagl^jJE" -1 


( 185 ) 


(186) 


where 


A j f°° S(oj)du ) 
7T JO OJ 2 T A j 


(187) 


So we’ve traded in a real matrix integral for a set of scalar, but often complex integrals. 
Here, we’ll carry out this program for cubic noise. 


The power spectrum in this case is (31), so for any given eigenvalue we have: 


* 



The 3 terms here may all be found in tables of indefinite integrals, when: 

/ \ "1 

• ( 4„ 2 + 3A 2 ) ton- 1 (j) - + w 2 + A 2 In (u> 2 + A 2 )^ (189) 


* = 


That this expression holds for complex A is readily checked by differentiation. On evalu- 
ating, this works out to: 


4> = 


2fl(0) /2w, 


h 


(¥) 


(190) 


where 


h{z) 


2 



tan ^ (z) — 2 



an even function of z. This looks very badly behaved near z — 0; but 
produced this Maclaurin expansion: 


(191) 

Matliematica 


h{z) =1- Jg 2 + w 2 - 


70 


105 




6 10 

1001 


+ ... 


(192) 


At the other end, the properties of the arctangent show that, as r — + oo, }j,{ z ) 
{7rsgn[3?(A)]}/2. The programs use the expansion for small z\ otherwise, standard for- 
mulas for complex logarithms and arctangents are employed. Note that if A is complex, 
its conjugate will also be present. Then, as <f>( A) is expressible as a power series in A, 
<f>(\) = 0(A); i.e., the function of the conjugate is equal to the conjugate of the function. 
This halves the arithmetic in such cases. All well and good; but the technique foundered 
on the numerical problems in obtaining complex eigenvectors. 
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D Appendix: Colored Noise Integral 


Consider the scalar integral 

2A uc f°° du ( 193 ) 

7T JO (u> 2 + A 2 )(w 2 + W 2 ) 

where A may be complex. If 9?(A) > 0, a contour integration along the real axis, and 
closed by the upper infinite half circle leads to 

J — (A I w c ) _1 O 94 ) 

while if 3?(A) < 0, a similar technique yields 

J = (A-Wc)" 1 (195) 

Both of these results have been checked by Mathematica. FinaHy, if 3ft(A) = 0, the integral 
is improper; but a limiting technique yields 

J = A( W ?-AV 0") 

Mathematica failed in this case, but [13], 3.264-1 provided verification. 


Now, for colored noise 


S{u>) = 2u; c «(0)(wM<) 


2\- 1 


and the integral (90) becomes 


$ 


2u> c R{0) 




7T 


w 2 + w 2 


Suppose the eigenvalue-eigenvector decomposition (183) is possible for Z. Then: 


2w c J?(0) 

7T 

2w c /J(0) 


A -V u; A 


f ( 

®B<iL g [a r 

7T L 


- .] E - 


(u> 2 + A 2 )(u> 2 + u; 2 ) 


= /t(0)Ediag|J(A)]E 


-l 


(197) 

(198) 


(199) 


In general, this is a mess. But if Z < 0: 

* = Jl(0)Ediag [(A - Wc )-»] E" 1 = H(0)E(A - t = Jl(0)(Z - ^I)" 1 (200) 


Similarly, if Z > 0, then 

$ = fi(0)(Z + w c I)- x (201) 

So, for any definite Z, the integration requires only a single matrix inversion, using real 
arithmetic. However, if Z is indefinite, the solution (199) may not be an improvement over 
direct numerical integration, that at least doesn’t require complex arithmetic. Fortunately, 
only Z < 0 can occur in the current filter theory. 
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